diff --git a/Code/Source/svFSI/ComMod.h b/Code/Source/svFSI/ComMod.h index 8ca9e7bb..f5bb6403 100644 --- a/Code/Source/svFSI/ComMod.h +++ b/Code/Source/svFSI/ComMod.h @@ -386,7 +386,7 @@ class stModelType /// @brief Fluid viscosity model type // -class viscModelType +class fluidViscModelType { public: @@ -409,6 +409,19 @@ class viscModelType double n = 0.0; }; +/// @brief Fluid viscosity model type +// +class solidViscModelType +{ + public: + + // Type of constitutive model for fluid viscosity + consts::SolidViscosityModelType viscType = consts::SolidViscosityModelType::viscType_NA; + + // Viscosity value + double mu = 0.0; +}; + /// @brief Domain type is to keep track with element belong to which domain /// and also different physical quantities // @@ -439,7 +452,10 @@ class dmnType stModelType stM; // Viscosity model for fluids - viscModelType visc; + fluidViscModelType fluid_visc; + + // Viscosity model for solids + solidViscModelType solid_visc; }; /// @brief Mesh adjacency (neighboring element for each element) diff --git a/Code/Source/svFSI/Parameters.cpp b/Code/Source/svFSI/Parameters.cpp index 72afc07d..bcea3dde 100644 --- a/Code/Source/svFSI/Parameters.cpp +++ b/Code/Source/svFSI/Parameters.cpp @@ -1069,43 +1069,43 @@ void VariableWallPropsParameters::set_values(tinyxml2::XMLElement* xml_elem) } ////////////////////////////////////////////////////////// -// ViscosityParameters // +// FluidViscosityParameters // ////////////////////////////////////////////////////////// -/// @brief Process parameters for various viscosity models. +/// @brief Process parameters for various fluid viscosity models. /// /// Define the XML element name for viscosiity parameters. -const std::string ViscosityParameters::xml_element_name_ = "Viscosity"; +const std::string FluidViscosityParameters::xml_element_name_ = "Viscosity"; -const std::string ViscosityParameters::CONSTANT_MODEL = "Constant"; -const std::string ViscosityParameters::CARREAU_YASUDA_MODEL = "Carreau-Yasuda"; -const std::string ViscosityParameters::CASSONS_MODEL = "Cassons"; +const std::string FluidViscosityParameters::CONSTANT_MODEL = "Constant"; +const std::string FluidViscosityParameters::CARREAU_YASUDA_MODEL = "Carreau-Yasuda"; +const std::string FluidViscosityParameters::CASSONS_MODEL = "Cassons"; -const std::set ViscosityParameters::model_names = { - ViscosityParameters::CONSTANT_MODEL, - ViscosityParameters::CARREAU_YASUDA_MODEL, - ViscosityParameters::CASSONS_MODEL +const std::set FluidViscosityParameters::model_names = { + FluidViscosityParameters::CONSTANT_MODEL, + FluidViscosityParameters::CARREAU_YASUDA_MODEL, + FluidViscosityParameters::CASSONS_MODEL }; -/// @brief Define a map to set parameters for each viscosity model. -using VpType = ViscosityParameters*; +/// @brief Define a map to set parameters for each fluid viscosity model. +using FVpType = FluidViscosityParameters*; using XmlType = tinyxml2::XMLElement*; -using SetViscosityParamMapType = std::map>; -SetViscosityParamMapType SetViscosityModelParamsMap = { - {ViscosityParameters::CARREAU_YASUDA_MODEL, [](VpType vp, XmlType params) -> void { vp->carreau_yasuda_model.set_values(params); }}, - {ViscosityParameters::CASSONS_MODEL, [](VpType vp, XmlType params) -> void { vp->cassons_model.set_values(params); }}, - {ViscosityParameters::CONSTANT_MODEL, [](VpType vp, XmlType params) -> void { vp->newtonian_model.set_values(params); }}, +using SetFluidViscosityParamMapType = std::map>; +SetFluidViscosityParamMapType SetFluidViscosityModelParamsMap = { + {FluidViscosityParameters::CARREAU_YASUDA_MODEL, [](FVpType vp, XmlType params) -> void { vp->carreau_yasuda_model.set_values(params); }}, + {FluidViscosityParameters::CASSONS_MODEL, [](FVpType vp, XmlType params) -> void { vp->cassons_model.set_values(params); }}, + {FluidViscosityParameters::CONSTANT_MODEL, [](FVpType vp, XmlType params) -> void { vp->newtonian_model.set_values(params); }}, }; -/// @brief Define a map to print parameters for each viscosity model. -using PrintViscosityParamaMapType = std::map>; -PrintViscosityParamaMapType PrintViscosityModelParamsMap = { - {ViscosityParameters::CARREAU_YASUDA_MODEL, [](VpType vp) -> void { vp->carreau_yasuda_model.print_parameters(); }}, - {ViscosityParameters::CASSONS_MODEL, [](VpType vp) -> void { vp->cassons_model.print_parameters(); }}, - {ViscosityParameters::CONSTANT_MODEL, [](VpType vp) -> void { vp->newtonian_model.print_parameters(); }}, +/// @brief Define a map to print parameters for each fluid viscosity model. +using PrintFluidViscosityParamaMapType = std::map>; +PrintFluidViscosityParamaMapType PrintFluidViscosityModelParamsMap = { + {FluidViscosityParameters::CARREAU_YASUDA_MODEL, [](FVpType vp) -> void { vp->carreau_yasuda_model.print_parameters(); }}, + {FluidViscosityParameters::CASSONS_MODEL, [](FVpType vp) -> void { vp->cassons_model.print_parameters(); }}, + {FluidViscosityParameters::CONSTANT_MODEL, [](FVpType vp) -> void { vp->newtonian_model.print_parameters(); }}, }; -ViscosityNewtonianParameters::ViscosityNewtonianParameters() +FluidViscosityNewtonianParameters::FluidViscosityNewtonianParameters() { // A parameter that must be defined. bool required = true; @@ -1113,24 +1113,24 @@ ViscosityNewtonianParameters::ViscosityNewtonianParameters() set_parameter("Value", 0.0, !required, constant_value); } -void ViscosityNewtonianParameters::print_parameters() +void FluidViscosityNewtonianParameters::print_parameters() { std::cout << constant_value.name_ << ": " << constant_value.value_ << std::endl; } -void ViscosityNewtonianParameters::set_values(tinyxml2::XMLElement* xml_elem) +void FluidViscosityNewtonianParameters::set_values(tinyxml2::XMLElement* xml_elem) { std::string error_msg = "Unknown Constitutive_model type=Newtonian XML element '"; using std::placeholders::_1; using std::placeholders::_2; std::function ftpr = - std::bind( &ViscosityNewtonianParameters::set_parameter_value, *this, _1, _2); + std::bind( &FluidViscosityNewtonianParameters::set_parameter_value, *this, _1, _2); xml_util_set_parameters(ftpr, xml_elem, error_msg); } -ViscosityCarreauYasudaParameters::ViscosityCarreauYasudaParameters() +FluidViscosityCarreauYasudaParameters::FluidViscosityCarreauYasudaParameters() { // A parameter that must be defined. bool required = true; @@ -1144,7 +1144,7 @@ ViscosityCarreauYasudaParameters::ViscosityCarreauYasudaParameters() set_parameter("Shear_rate_tensor_exponent", 0.0, !required, shear_rate_tensor_exponent); } -void ViscosityCarreauYasudaParameters::print_parameters() +void FluidViscosityCarreauYasudaParameters::print_parameters() { std::cout << limiting_high_shear_rate_viscosity.name_ << ": " << limiting_high_shear_rate_viscosity.value_ << std::endl; std::cout << limiting_low_shear_rate_viscosity.name_ << ": " << limiting_low_shear_rate_viscosity.value_ << std::endl; @@ -1153,18 +1153,18 @@ void ViscosityCarreauYasudaParameters::print_parameters() std::cout << shear_rate_tensor_multipler.name_ << ": " << shear_rate_tensor_multipler.value_ << std::endl; } -void ViscosityCarreauYasudaParameters::set_values(tinyxml2::XMLElement* xml_elem) +void FluidViscosityCarreauYasudaParameters::set_values(tinyxml2::XMLElement* xml_elem) { std::string error_msg = "Unknown Constitutive_model type=CarreauYasuda XML element '"; using std::placeholders::_1; using std::placeholders::_2; std::function ftpr = - std::bind( &ViscosityCarreauYasudaParameters::set_parameter_value, *this, _1, _2); + std::bind( &FluidViscosityCarreauYasudaParameters::set_parameter_value, *this, _1, _2); xml_util_set_parameters(ftpr, xml_elem, error_msg); } -ViscosityCassonsParameters::ViscosityCassonsParameters() +FluidViscosityCassonsParameters::FluidViscosityCassonsParameters() { // A parameter that must be defined. bool required = true; @@ -1174,26 +1174,26 @@ ViscosityCassonsParameters::ViscosityCassonsParameters() set_parameter("Yield_stress_parameter", 0.0, !required, yield_stress); } -void ViscosityCassonsParameters::print_parameters() +void FluidViscosityCassonsParameters::print_parameters() { std::cout << asymptotic_viscosity.name_ << ": " << asymptotic_viscosity.value_ << std::endl; std::cout << low_shear_rate_threshold.name_ << ": " << low_shear_rate_threshold.value_ << std::endl; std::cout << yield_stress.name_ << ": " << yield_stress.value_ << std::endl; } -void ViscosityCassonsParameters::set_values(tinyxml2::XMLElement* xml_elem) +void FluidViscosityCassonsParameters::set_values(tinyxml2::XMLElement* xml_elem) { std::string error_msg = "Unknown Constitutive_model type=Cassons XML element '"; using std::placeholders::_1; using std::placeholders::_2; std::function ftpr = - std::bind( &ViscosityCassonsParameters::set_parameter_value, *this, _1, _2); + std::bind( &FluidViscosityCassonsParameters::set_parameter_value, *this, _1, _2); xml_util_set_parameters(ftpr, xml_elem, error_msg); } -ViscosityParameters::ViscosityParameters() +FluidViscosityParameters::FluidViscosityParameters() { // A parameter that must be defined. bool required = true; @@ -1202,7 +1202,7 @@ ViscosityParameters::ViscosityParameters() model = Parameter("model", "", required); } -void ViscosityParameters::print_parameters() +void FluidViscosityParameters::print_parameters() { std::cout << std::endl; std::cout << "--------------------" << std::endl; @@ -1210,11 +1210,11 @@ void ViscosityParameters::print_parameters() std::cout << "--------------------" << std::endl; std::cout << model.name() << ": '" << model.value() << "'" << std::endl; - // Print parameters for the given viscosity model. - PrintViscosityModelParamsMap[model.value_](this); + // Print parameters for the given fluid_viscosity model. + PrintFluidViscosityModelParamsMap[model.value_](this); } -void ViscosityParameters::set_values(tinyxml2::XMLElement* xml_elem) +void FluidViscosityParameters::set_values(tinyxml2::XMLElement* xml_elem) { using namespace tinyxml2; @@ -1226,14 +1226,139 @@ void ViscosityParameters::set_values(tinyxml2::XMLElement* xml_elem) } model.set(std::string(smodel)); - // Check viscosity model name. + // Check fluid_viscosity model name. if (model_names.count(model.value()) == 0) { - throw std::runtime_error("Unknown viscosity model '" + model.value() + - " in '" + xml_elem->Name() + "'."); + throw std::runtime_error("Unknown fluid viscosity model '" + model.value() + + "' in '" + xml_elem->Name() + "'."); } - // Set parameters for the given viscosity model. - SetViscosityModelParamsMap[model.value()](this, xml_elem); + // Set parameters for the given fluid_viscosity model. + SetFluidViscosityModelParamsMap[model.value()](this, xml_elem); +} + +////////////////////////////////////////////////////////// +// SolidViscosityParameters // +////////////////////////////////////////////////////////// + +/// @brief Process parameters for various solid viscosity models. +/// +/// Define the XML element name for viscosiity parameters. +const std::string SolidViscosityParameters::xml_element_name_ = "Viscosity"; + +const std::string SolidViscosityParameters::NEWTONIAN_MODEL = "Newtonian"; +const std::string SolidViscosityParameters::POTENTIAL_MODEL = "Potential"; + +const std::set SolidViscosityParameters::model_names = { + SolidViscosityParameters::NEWTONIAN_MODEL, + SolidViscosityParameters::POTENTIAL_MODEL +}; + +/// @brief Define a map to set parameters for each solid viscosity model. +using SVpType = SolidViscosityParameters*; +using XmlType = tinyxml2::XMLElement*; +using SetSolidViscosityParamMapType = std::map>; +SetSolidViscosityParamMapType SetSolidViscosityModelParamsMap = { + {SolidViscosityParameters::NEWTONIAN_MODEL, [](SVpType vp, XmlType params) -> void { vp->newtonian_model.set_values(params); }}, + {SolidViscosityParameters::POTENTIAL_MODEL, [](SVpType vp, XmlType params) -> void { vp->potential_model.set_values(params); }}, +}; + +/// @brief Define a map to print parameters for each solid viscosity model. +using PrintSolidViscosityParamaMapType = std::map>; +PrintSolidViscosityParamaMapType PrintSolidViscosityModelParamsMap = { + {SolidViscosityParameters::NEWTONIAN_MODEL, [](SVpType vp) -> void { vp->newtonian_model.print_parameters(); }}, + {SolidViscosityParameters::POTENTIAL_MODEL, [](SVpType vp) -> void { vp->potential_model.print_parameters(); }}, +}; + +SolidViscosityNewtonianParameters::SolidViscosityNewtonianParameters() +{ + // A parameter that must be defined. + bool required = true; + + set_parameter("Value", 0.0, !required, constant_value); +} + +void SolidViscosityNewtonianParameters::print_parameters() +{ + std::cout << constant_value.name_ << ": " << constant_value.value_ << std::endl; +} + +void SolidViscosityNewtonianParameters::set_values(tinyxml2::XMLElement* xml_elem) +{ + std::string error_msg = "Unknown Constitutive_model type=Newtonian XML element '"; + + using std::placeholders::_1; + using std::placeholders::_2; + std::function ftpr = + std::bind( &SolidViscosityNewtonianParameters::set_parameter_value, *this, _1, _2); + + xml_util_set_parameters(ftpr, xml_elem, error_msg); +} + +SolidViscosityPotentialParameters::SolidViscosityPotentialParameters() +{ + // A parameter that must be defined. + bool required = true; + + set_parameter("Value", 0.0, !required, constant_value); +} + +void SolidViscosityPotentialParameters::print_parameters() +{ + std::cout << constant_value.name_ << ": " << constant_value.value_ << std::endl; +} + +void SolidViscosityPotentialParameters::set_values(tinyxml2::XMLElement* xml_elem) +{ + std::string error_msg = "Unknown Constitutive_model type=Potential XML element '"; + using std::placeholders::_1; + using std::placeholders::_2; + std::function ftpr = + std::bind( &SolidViscosityPotentialParameters::set_parameter_value, *this, _1, _2); + + xml_util_set_parameters(ftpr, xml_elem, error_msg); +} + +SolidViscosityParameters::SolidViscosityParameters() +{ + // A parameter that must be defined. + bool required = true; + + // Stores model from the XML element + model = Parameter("model", "", required); +} + +void SolidViscosityParameters::print_parameters() +{ + std::cout << std::endl; + std::cout << "--------------------" << std::endl; + std::cout << "Viscosity Parameters" << std::endl; + std::cout << "--------------------" << std::endl; + std::cout << model.name() << ": '" << model.value() << "'" << std::endl; + + // Print parameters for the given solid_viscosity model. + PrintSolidViscosityModelParamsMap[model.value_](this); +} + +void SolidViscosityParameters::set_values(tinyxml2::XMLElement* xml_elem) +{ + using namespace tinyxml2; + + const char* smodel; + auto result = xml_elem->QueryStringAttribute("model", &smodel); + + if (smodel == nullptr) { + throw std::runtime_error("No MODEL given in the XML element."); + } + model.set(std::string(smodel)); + + // Check solid viscosity model name. + if (model_names.count(model.value()) == 0) { + throw std::runtime_error("Unknown solid viscosity model '" + model.value() + + "' in '" + xml_elem->Name() + "'."); + } + + // Set parameters for the given solid viscosity model. + SetSolidViscosityModelParamsMap[model.value()](this, xml_elem); } ////////////////////////////////////////////////////////// @@ -1299,7 +1424,6 @@ DomainParameters::DomainParameters() set_parameter("Relative_tolerance", 1e-4, !required, relative_tolerance); set_parameter("Shell_thickness", 0.0, !required, shell_thickness); set_parameter("Solid_density", 0.5, !required, solid_density); - set_parameter("Solid_viscosity", 0.0, !required, solid_viscosity); set_parameter("Source_term", 0.0, !required, source_term); set_parameter("Time_step_for_integration", 0.0, !required, time_step_for_integration); } @@ -1323,7 +1447,9 @@ void DomainParameters::print_parameters() stimulus.print_parameters(); - viscosity.print_parameters(); + fluid_viscosity.print_parameters(); + + solid_viscosity.print_parameters(); } void DomainParameters::set_values(tinyxml2::XMLElement* domain_elem) @@ -1354,8 +1480,16 @@ void DomainParameters::set_values(tinyxml2::XMLElement* domain_elem) } else if (name == StimulusParameters::xml_element_name_) { stimulus.set_values(item); - } else if (name == ViscosityParameters::xml_element_name_) { - viscosity.set_values(item); + } else if (name == FluidViscosityParameters::xml_element_name_ || name == SolidViscosityParameters::xml_element_name_) { + auto eq_type = consts::equation_name_to_type.at(equation.value()); + if (eq_type == consts::EquationType::phys_fluid || eq_type == consts::EquationType::phys_CMM || eq_type == consts::EquationType::phys_stokes) { + fluid_viscosity.set_values(item); + } else if (eq_type == consts::EquationType::phys_struct || eq_type == consts::EquationType::phys_ustruct) { + solid_viscosity.set_values(item); + } + else { + throw std::runtime_error("Viscosity model not supported for equation '" + equation.value() + "'."); + } } else if (item->GetText() != nullptr) { auto value = item->GetText(); @@ -1768,8 +1902,15 @@ void EquationParameters::set_values(tinyxml2::XMLElement* eq_elem) } else if (name == StimulusParameters::xml_element_name_) { default_domain->stimulus.set_values(item); - } else if (name == ViscosityParameters::xml_element_name_) { - default_domain->viscosity.set_values(item); + } else if (name == FluidViscosityParameters::xml_element_name_ || name == SolidViscosityParameters::xml_element_name_) { + auto eq_type = consts::equation_name_to_type.at(type.value()); + if (eq_type == consts::EquationType::phys_fluid || eq_type == consts::EquationType::phys_CMM || eq_type == consts::EquationType::phys_stokes) { + default_domain->fluid_viscosity.set_values(item); + } else if (eq_type == consts::EquationType::phys_struct || eq_type == consts::EquationType::phys_ustruct) { + default_domain->solid_viscosity.set_values(item); + } else { + throw std::runtime_error("Viscosity model not supported for equation '" + type.value() + "'."); + } } else if (name == ECGLeadsParameters::xml_element_name_) { ecg_leads.set_values(item); diff --git a/Code/Source/svFSI/Parameters.h b/Code/Source/svFSI/Parameters.h index 25f1d79b..4e1a7830 100644 --- a/Code/Source/svFSI/Parameters.h +++ b/Code/Source/svFSI/Parameters.h @@ -791,25 +791,25 @@ class VariableWallPropsParameters : public ParameterLists ////////////////////////////////////////////////////////// -// Viscosity // +// FluidViscosity // ////////////////////////////////////////////////////////// // The following classes are used to store parameters for -// various viscosity models. +// various fluid viscosity models. -class ViscosityNewtonianParameters : public ParameterLists +class FluidViscosityNewtonianParameters : public ParameterLists { public: - ViscosityNewtonianParameters(); + FluidViscosityNewtonianParameters(); void print_parameters(); void set_values(tinyxml2::XMLElement* equation_params); Parameter constant_value; }; -class ViscosityCarreauYasudaParameters : public ParameterLists +class FluidViscosityCarreauYasudaParameters : public ParameterLists { public: - ViscosityCarreauYasudaParameters(); + FluidViscosityCarreauYasudaParameters(); void print_parameters(); void set_values(tinyxml2::XMLElement* xml_elem); @@ -820,10 +820,10 @@ class ViscosityCarreauYasudaParameters : public ParameterLists Parameter shear_rate_tensor_exponent; }; -class ViscosityCassonsParameters : public ParameterLists +class FluidViscosityCassonsParameters : public ParameterLists { public: - ViscosityCassonsParameters(); + FluidViscosityCassonsParameters(); void print_parameters(); void set_values(tinyxml2::XMLElement* xml_elem); Parameter asymptotic_viscosity; @@ -831,10 +831,10 @@ class ViscosityCassonsParameters : public ParameterLists Parameter low_shear_rate_threshold; }; -class ViscosityParameters : public ParameterLists +class FluidViscosityParameters : public ParameterLists { public: - ViscosityParameters(); + FluidViscosityParameters(); static const std::string xml_element_name_; @@ -848,11 +848,57 @@ class ViscosityParameters : public ParameterLists Parameter model; - ViscosityNewtonianParameters newtonian_model; - ViscosityCarreauYasudaParameters carreau_yasuda_model; - ViscosityCassonsParameters cassons_model; + FluidViscosityNewtonianParameters newtonian_model; + FluidViscosityCarreauYasudaParameters carreau_yasuda_model; + FluidViscosityCassonsParameters cassons_model; }; +////////////////////////////////////////////////////////// +// SolidViscosity // +////////////////////////////////////////////////////////// + +// The following classes are used to store parameters for +// various solid viscosity models. + +class SolidViscosityNewtonianParameters : public ParameterLists +{ + public: + SolidViscosityNewtonianParameters(); + void print_parameters(); + void set_values(tinyxml2::XMLElement* equation_params); + Parameter constant_value; +}; + +class SolidViscosityPotentialParameters : public ParameterLists +{ + public: + SolidViscosityPotentialParameters(); + void print_parameters(); + void set_values(tinyxml2::XMLElement* equation_params); + Parameter constant_value; +}; + +class SolidViscosityParameters : public ParameterLists +{ + public: + SolidViscosityParameters(); + + static const std::string xml_element_name_; + + static const std::string NEWTONIAN_MODEL; + static const std::string POTENTIAL_MODEL; + static const std::set model_names; + + void print_parameters(); + void set_values(tinyxml2::XMLElement* xml_elem); + + Parameter model; + + SolidViscosityNewtonianParameters newtonian_model; + SolidViscosityPotentialParameters potential_model; +}; + + /// @brief The LinearAlgebraParameters class stores parameters for /// the 'Linear_algebra' XML element. class LinearAlgebraParameters : public ParameterLists @@ -1011,7 +1057,8 @@ class DomainParameters : public ParameterLists ConstitutiveModelParameters constitutive_model; FiberReinforcementStressParameters fiber_reinforcement_stress; StimulusParameters stimulus; - ViscosityParameters viscosity; + FluidViscosityParameters fluid_viscosity; + SolidViscosityParameters solid_viscosity; // Attributes. Parameter id; @@ -1060,7 +1107,6 @@ class DomainParameters : public ParameterLists Parameter shell_thickness; Parameter solid_density; - Parameter solid_viscosity; Parameter source_term; Parameter time_step_for_integration; }; @@ -1199,7 +1245,9 @@ class EquationParameters : public ParameterLists VariableWallPropsParameters variable_wall_properties; - ViscosityParameters viscosity; + FluidViscosityParameters fluid_viscosity; + + SolidViscosityParameters solid_viscosity; ECGLeadsParameters ecg_leads; }; diff --git a/Code/Source/svFSI/README.md b/Code/Source/svFSI/README.md index c9f84a50..e1b4b622 100644 --- a/Code/Source/svFSI/README.md +++ b/Code/Source/svFSI/README.md @@ -1048,7 +1048,7 @@ Distribute boundary condition data to all processors. -

dist_visc_model(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, viscModelType& lVis)

+

dist_visc_model(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, fluidViscModelType& lVis)

[distribute.cpp](https://github.com/ktbolt/svFSI/blob/Implement-svFSI-using-cpp_19/Code/Source/svFSI_cinterface/distribute.cpp) diff --git a/Code/Source/svFSI/consts.cpp b/Code/Source/svFSI/consts.cpp index 0ee1f2cf..bd2e521c 100644 --- a/Code/Source/svFSI/consts.cpp +++ b/Code/Source/svFSI/consts.cpp @@ -126,6 +126,21 @@ const std::map fluid_viscosity_model_name_t }; +/// @brief Map for solid viscosity model string to SolidViscosityModelType. +// +const std::map solid_viscosity_model_name_to_type +{ + {"Newtonian", SolidViscosityModelType::viscType_Newtonian}, + {"newtonian", SolidViscosityModelType::viscType_Newtonian}, + {"Newt", SolidViscosityModelType::viscType_Newtonian}, + {"newt", SolidViscosityModelType::viscType_Newtonian}, + + {"Potential", SolidViscosityModelType::viscType_Potential}, + {"potential", SolidViscosityModelType::viscType_Potential}, + {"Pot", SolidViscosityModelType::viscType_Potential}, + {"pot", SolidViscosityModelType::viscType_Potential}, +}; + /// @brief Map number of element nodes to element type. /// diff --git a/Code/Source/svFSI/consts.h b/Code/Source/svFSI/consts.h index 8d7e9720..f2790ba3 100644 --- a/Code/Source/svFSI/consts.h +++ b/Code/Source/svFSI/consts.h @@ -398,19 +398,18 @@ enum class PhysicalProperyType NA = 0, fluid_density = 1, solid_density = 2, - solid_viscosity = 3, - elasticity_modulus = 4, - poisson_ratio = 5, - conductivity = 6, - f_x = 7, // internal force x - f_y = 8, // internal force y - f_z = 9, // internal force z - backflow_stab = 10, // stabilization coeff. for backflow divergence - source_term = 11, // external source - damping = 12, - shell_thickness = 13, - ctau_M = 14, // stabilization coeffs. for USTRUCT (momentum, continuity) - ctau_C = 15 + elasticity_modulus = 3, + poisson_ratio = 4, + conductivity = 5, + f_x = 6, // internal force x + f_y = 7, // internal force y + f_z = 8, // internal force z + backflow_stab = 9, // stabilization coeff. for backflow divergence + source_term = 10, // external source + damping = 11, + shell_thickness = 12, + ctau_M = 13, // stabilization coeffs. for USTRUCT (momentum, continuity) + ctau_C = 14 }; enum class PreconditionerType @@ -460,6 +459,16 @@ enum class FluidViscosityModelType /// Map for fluid viscosity model string to FluidViscosityModelType. extern const std::map fluid_viscosity_model_name_to_type; +enum class SolidViscosityModelType +{ + viscType_NA = 695, + viscType_Newtonian = 694, + viscType_Potential = 693 +}; + +/// Map for solid viscosity model string to SolidViscosityModelType. +extern const std::map solid_viscosity_model_name_to_type; + /// Template for printing enum class types as an int. template std::ostream& operator<<(typename std::enable_if::value, std::ostream>::type& stream, const T& e) diff --git a/Code/Source/svFSI/distribute.cpp b/Code/Source/svFSI/distribute.cpp index d67f0f84..93532359 100644 --- a/Code/Source/svFSI/distribute.cpp +++ b/Code/Source/svFSI/distribute.cpp @@ -1061,12 +1061,13 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std:: if ((dmn.phys == EquationType::phys_struct) || (dmn.phys == EquationType::phys_ustruct)) { dist_mat_consts(com_mod, cm_mod, cm, dmn.stM); + dist_solid_visc_model(com_mod, cm_mod, cm, dmn.solid_visc); } if ((dmn.phys == EquationType::phys_fluid) || (dmn.phys == EquationType::phys_stokes) || (dmn.phys == EquationType::phys_CMM && !com_mod.cmmInit)) { - dist_visc_model(com_mod, cm_mod, cm, dmn.visc); + dist_fluid_visc_model(com_mod, cm_mod, cm, dmn.fluid_visc); } } @@ -1224,7 +1225,7 @@ void dist_mat_consts(const ComMod& com_mod, const CmMod& cm_mod, const cmType& c } -void dist_visc_model(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, viscModelType& lVis) +void dist_fluid_visc_model(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, fluidViscModelType& lVis) { using namespace consts; @@ -1236,6 +1237,14 @@ void dist_visc_model(const ComMod& com_mod, const CmMod& cm_mod, const cmType& c cm.bcast(cm_mod, &lVis.n); } +void dist_solid_visc_model(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, solidViscModelType& lVis) +{ + using namespace consts; + + cm.bcast_enum(cm_mod, &lVis.viscType); + cm.bcast(cm_mod, &lVis.mu); +} + void part_face(Simulation* simulation, mshType& lM, faceType& lFa, faceType& gFa, Vector& gmtl) { diff --git a/Code/Source/svFSI/distribute.h b/Code/Source/svFSI/distribute.h index bd831f79..33cfabd2 100644 --- a/Code/Source/svFSI/distribute.h +++ b/Code/Source/svFSI/distribute.h @@ -45,7 +45,9 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std:: void dist_mat_consts(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, stModelType& lStM); -void dist_visc_model(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, viscModelType& lVis); +void dist_fluid_visc_model(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, fluidViscModelType& lVis); + +void dist_solid_visc_model(const ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, solidViscModelType& lVis); void part_face(Simulation* simulation, mshType& lM, faceType& lFa, faceType& gFa, Vector& gmtl); diff --git a/Code/Source/svFSI/fluid.cpp b/Code/Source/svFSI/fluid.cpp index 2318e665..c9e35109 100644 --- a/Code/Source/svFSI/fluid.cpp +++ b/Code/Source/svFSI/fluid.cpp @@ -1950,20 +1950,20 @@ void get_viscosity(const ComMod& com_mod, const dmnType& lDmn, double& gamma, do double mu_i, mu_o, lam, a, n, T1, T2; - switch (lDmn.visc.viscType) { + switch (lDmn.fluid_visc.viscType) { case FluidViscosityModelType::viscType_Const: - mu = lDmn.visc.mu_i; + mu = lDmn.fluid_visc.mu_i; mu_s = mu; mu_x = 0.0; break; case FluidViscosityModelType::viscType_CY: - mu_i = lDmn.visc.mu_i; - mu_o = lDmn.visc.mu_o; - lam = lDmn.visc.lam; - a = lDmn.visc.a; - n = lDmn.visc.n; + mu_i = lDmn.fluid_visc.mu_i; + mu_o = lDmn.fluid_visc.mu_o; + lam = lDmn.fluid_visc.lam; + a = lDmn.fluid_visc.a; + n = lDmn.fluid_visc.n; T1 = 1.0 + pow(lam*gamma, a); T2 = pow(T1,((n-1.0)/a)); @@ -1976,9 +1976,9 @@ void get_viscosity(const ComMod& com_mod, const dmnType& lDmn, double& gamma, do break; case FluidViscosityModelType::viscType_Cass: - mu_i = lDmn.visc.mu_i; - mu_o = lDmn.visc.mu_o; - lam = lDmn.visc.lam; + mu_i = lDmn.fluid_visc.mu_i; + mu_o = lDmn.fluid_visc.mu_o; + lam = lDmn.fluid_visc.lam; if (gamma < lam) { mu_o = mu_o / sqrt(lam); diff --git a/Code/Source/svFSI/mat_models_carray.h b/Code/Source/svFSI/mat_models_carray.h index 8946fcd6..85f80d31 100644 --- a/Code/Source/svFSI/mat_models_carray.h +++ b/Code/Source/svFSI/mat_models_carray.h @@ -39,6 +39,7 @@ #include "Tensor4.h" #include "mat_fun.h" +#include "mat_fun_carray.h" namespace mat_models_carray { @@ -1134,6 +1135,278 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn cc_to_voigt_carray(CC, Dm); } + +/** + * @brief Get the viscous PK2 stress and corresponding tangent matrix contributions for a solid + * with a viscous pseudo-potential model. + * This is defined by a viscous pseuo-potential + * Psi = mu/2 * tr(E_dot^2) + * The viscous 2nd Piola-Kirchhoff stress is given by + * Svis = dPsi/dE_dot + * = mu * E_dot + * = mu * 1/2 * F^T * (grad(v) + grad(v)^T) * F + * = mu * 1/2 * ( (F^T * Grad(v)) + (F^T * Grad(v))^T ) + * + * @tparam nsd Number of spatial dimensions + * @param mu Solid viscosity parameter + * @param eNoN Number of nodes in an element + * @param Nx Shape function gradient w.r.t. reference configuration coordinates (dN/dX) + * @param vx Velocity gradient matrix w.r.t reference configuration coordinates (dv/dX) + * @param F Deformation gradient matrix + * @param Svis Viscous 2nd Piola-Kirchhoff stress matrix + * @param Kvis_u Viscous tangent matrix contribution due to displacement + * @param Kvis_v Visous tangent matrix contribution due to velocity + */ +template +void get_visc_stress_pot(const double mu, const int eNoN, const Array& Nx, const double vx[nsd][nsd], const double F[nsd][nsd], + Array& Svis, Array3& Kvis_u, Array3& Kvis_v) { + + + // Initialize Svis, Kvis_u, Kvis_v to zero + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < nsd; j++) { + Svis(i,j) = 0.0; + for (int a = 0; a < eNoN; a++) { + for (int b = 0; b < eNoN; b++) { + Kvis_u(i*nsd+j,a,b) = 0.0; + Kvis_v(i*nsd+j,a,b) = 0.0; + } + } + } + } + + + // Required intermediate terms for stress and tangent + double Ft[nsd][nsd] = {0}, vxt[nsd][nsd] = {0}, F_Ft[nsd][nsd] = {0}, Ft_vx[nsd][nsd] = {0}, F_vxt[nsd][nsd] = {0}; + mat_fun_carray::transpose(F, Ft); + mat_fun_carray::mat_mul(F, Ft, F_Ft); + mat_fun_carray::mat_mul(Ft, vx, Ft_vx); + mat_fun_carray::transpose(vx, vxt); + mat_fun_carray::mat_mul(F, vxt, F_vxt); + + //double F_Nx[nsd][eNoN] = {0}, vx_Nx[nsd][eNoN] = {0}; + Array F_Nx(nsd,eNoN), vx_Nx(nsd,eNoN); + + for (int a = 0; a < eNoN; ++a) { + for (int i = 0; i < nsd; ++i) { + for (int j = 0; j < nsd; ++j) { + F_Nx(i,a) += F[i][j] * Nx(j,a); + vx_Nx(i,a) += vx[i][j] * Nx(j,a); + } + } + } + + // 2nd Piola-Kirchhoff stress due to viscosity + // Svis = mu * 1/2 * ( (F^T * dv/dX) + (F^T * dv/dX)^T ) + double Ft_vx_symm[nsd][nsd] = {0}; + mat_fun_carray::mat_symm(Ft_vx, Ft_vx_symm); + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < nsd; j++) { + Svis(i,j) = mu * Ft_vx_symm[i][j]; + } + } + + // Tangent matrix contributions due to viscosity + for (int b = 0; b < eNoN; ++b) { + for (int a = 0; a < eNoN; ++a) { + double Nx_Nx = 0.0; + for (int i = 0; i < nsd; ++i) { + Nx_Nx += Nx(i,a) * Nx(i,b); + } + + for (int i = 0; i < nsd; ++i) { + for (int j = 0; j < nsd; ++j) { + int ii = i * nsd + j; + Kvis_u(ii,a,b) = 0.5 * mu * (F_Nx(i,b) * vx_Nx(j,a) + Nx_Nx * F_vxt[i][j]); + Kvis_v(ii,a,b) = 0.5 * mu * (Nx_Nx * F_Ft[i][j] + F_Nx(i,b) * F_Nx(j,a)); + } + } + } + } +} + +/** + * @brief Get the viscous PK2 stress and corresponding tangent matrix contributions for a solid + * with a Newtonian fluid-like viscosity model. + * The viscous deviatoric Cauchy stress is given by + * sigma_vis_dev = 2 * mu * d_dev + * where d_dev = 1/2 * (grad(v) + grad(v)^T) - 1/3 * (div(v)) * I + * The viscous 2nd Piola-Kirchhoff stress is given by a pull-back operation + * Svis = 2 * mu * J * F^-1 * d_dev * F^-T + * + * Note, there is likely an error/bug in the tangent contributions that leads to suboptimal nonlinear convergence + * + * @tparam nsd Number of spatial dimensions + * @param mu Solid viscosity parameter + * @param eNoN Number of nodes in an element + * @param Nx Shape function gradient w.r.t. reference configuration coordinates (dN/dX) + * @param vx Velocity gradient matrix w.r.t reference configuration coordinates (dv/dX) + * @param F Deformation gradient matrix + * @param Svis Viscous 2nd Piola-Kirchhoff stress matrix + * @param Kvis_u Viscous tangent matrix contribution due to displacement + * @param Kvis_v Visous tangent matrix contribution due to velocity + */ +template +void get_visc_stress_newt(const double mu, const int eNoN, const Array& Nx, const double vx[nsd][nsd], const double F[nsd][nsd], + Array& Svis, Array3& Kvis_u, Array3& Kvis_v) { + + + // Initialize Svis, Kvis_u, Kvis_v to zero + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < nsd; j++) { + Svis(i,j) = 0.0; + for (int a = 0; a < eNoN; a++) { + for (int b = 0; b < eNoN; b++) { + Kvis_u(i*nsd+j,a,b) = 0.0; + Kvis_v(i*nsd+j,a,b) = 0.0; + } + } + } + } + + // Get identity matrix, Jacobian, and F^-1 + double Idm[nsd][nsd] = {0}; + mat_fun_carray::mat_id(Idm); + double J = mat_fun_carray::mat_det(F); + double Fi[nsd][nsd] = {0}; + mat_fun_carray::mat_inv(F, Fi); + + // Required intermediate terms for stress and tangent + double vx_Fi[nsd][nsd] = {0}, vx_Fi_symm[nsd][nsd] = {0}, ddev[nsd][nsd] = {0}; + // vx_Fi: Velocity gradient in current configuration + mat_fun_carray::mat_mul(vx, Fi, vx_Fi); + mat_fun_carray::mat_symm(vx_Fi, vx_Fi_symm); + // ddev: Deviatoric part of rate of strain tensor + mat_fun_carray::mat_dev(vx_Fi_symm, ddev); + //double Nx_Fi[nsd][eNoN] = {0}, ddev_Nx_Fi[nsd][eNoN] = {0}, vx_Fi_Nx_Fi[nsd][eNoN] = {0}; + Array Nx_Fi(nsd,eNoN), ddev_Nx_Fi(nsd,eNoN), vx_Fi_Nx_Fi(nsd,eNoN); + for (int a = 0; a < eNoN; ++a) { + for (int i = 0; i < nsd; ++i) { + for (int j = 0; j < nsd; ++j) { + Nx_Fi(i,a) += Nx(j,a) * Fi[j][i]; + } + } + for (int i = 0; i < nsd; ++i) { + for (int j = 0; j < nsd; ++j) { + ddev_Nx_Fi(i,a) += ddev[i][j] * Nx_Fi(j,a); + vx_Fi_Nx_Fi(i,a) += vx_Fi[j][i] * Nx_Fi(j,a); + } + } + } + + // 2nd Piola-Kirchhoff stress due to viscosity + // Svis = 2 * mu * J * F^-1 * d_dev * F^-T + double Fit[nsd][nsd] = {0}; + mat_fun_carray::transpose(Fi, Fit); + double ddev_Fit[nsd][nsd] = {0}; + mat_fun_carray::mat_mul(ddev, Fit, ddev_Fit); + double Fi_ddev_Fit[nsd][nsd] = {0}; + mat_fun_carray::mat_mul(Fi, ddev_Fit, Fi_ddev_Fit); + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < nsd; j++) { + Svis(i,j) = 2.0 * mu * J * Fi_ddev_Fit[i][j]; + } + } + + // Tangent matrix contributions due to viscosity + double r2d = 2.0 / nsd; + for (int b = 0; b < eNoN; ++b) { + for (int a = 0; a < eNoN; ++a) { + double Nx_Fi_Nx_Fi = 0.0; + for (int i = 0; i < nsd; ++i) { + Nx_Fi_Nx_Fi += Nx_Fi(i,a) * Nx_Fi(i,b); + } + + for (int i = 0; i < nsd; ++i) { + for (int j = 0; j < nsd; ++j) { + int ii = i * nsd + j; + + // Derivative of the residual w.r.t displacement + Kvis_u(ii,a,b) = mu * J * (2.0 * + (ddev_Nx_Fi(i,a) * Nx_Fi(j,b) - ddev_Nx_Fi(i,b) * Nx_Fi(j,a)) - + (Nx_Fi_Nx_Fi * vx_Fi[i][j] + Nx_Fi(i,b) * vx_Fi_Nx_Fi(j,a) - + r2d * Nx_Fi(i,a) * vx_Fi_Nx_Fi(j,b))); + + // Derivative of the residual w.r.t velocity + Kvis_v(ii,a,b) = mu * J * (Nx_Fi_Nx_Fi * Idm[i][j] + + Nx_Fi(i,b) * Nx_Fi(j,a) - r2d * Nx_Fi(i,a) * Nx_Fi(j,b)); + } + } + } + } +} + + +/** + * @brief Get the solid viscous PK2 stress and corresponding tangent matrix contributions + * Calls the appropriate function based on the viscosity type, either viscous + * pseudo-potential or Newtonian viscosity model. + * + * @tparam nsd Number of spatial dimensions + * @param[in] lDmn Domain object + * @param[in] eNoN Number of nodes in an element + * @param[in] Nx Shape function gradient w.r.t. reference configuration coordinates (dN/dX) + * @param[in] vx Velocity gradient matrix w.r.t reference configuration coordinates (dv/dX) + * @param[in] F Deformation gradient matrix + * @param[out] Svis Viscous 2nd Piola-Kirchhoff stress matrix + * @param[out] Kvis_u Viscous tangent matrix contribution due to displacement + * @param[out] Kvis_v Viscous tangent matrix contribution due to velocity + */ +template +void get_visc_stress_and_tangent(const dmnType& lDmn, const int eNoN, const Array& Nx, const double vx[nsd][nsd], const double F[nsd][nsd], + Array& Svis, Array3& Kvis_u, Array3& Kvis_v) { + + switch (lDmn.solid_visc.viscType) { + case consts::SolidViscosityModelType::viscType_Newtonian: + get_visc_stress_newt(lDmn.solid_visc.mu, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v); + break; + + case consts::SolidViscosityModelType::viscType_Potential: + get_visc_stress_pot(lDmn.solid_visc.mu, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v); + break; + } +} + +/** + * @brief Get the solid viscous PK2 stress and corresponding tangent matrix contributions + * Calls the appropriate function based on the viscosity type, either viscous + * pseudo-potential or Newtonian viscosity model. + * + * Same as above, except takes vx and F as Array objects instead of C-style arrays. + * + * @tparam nsd Number of spatial dimensions + * @param[in] lDmn Domain object + * @param[in] eNoN Number of nodes in an element + * @param[in] Nx Shape function gradient w.r.t. reference configuration coordinates (dN/dX) + * @param[in] vx Velocity gradient matrix w.r.t reference configuration coordinates (dv/dX) + * @param[in] F Deformation gradient matrix + * @param[out] Svis Viscous 2nd Piola-Kirchhoff stress matrix + * @param[out] Kvis_u Viscous tangent matrix contribution due to displacement + * @param[out] Kvis_v Viscous tangent matrix contribution due to velocity + */ +template +void get_visc_stress_and_tangent(const dmnType& lDmn, const int eNoN, const Array& Nx, Array& vx_Array, Array& F_Array, + Array& Svis, Array3& Kvis_u, Array3& Kvis_v) { + + // Convert vx_Array and F_Array to C-style arrays + double vx[nsd][nsd], F[nsd][nsd]; + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < nsd; j++) { + vx[i][j] = vx_Array(i,j); + F[i][j] = F_Array(i,j); + } + } + switch (lDmn.solid_visc.viscType) { + case consts::SolidViscosityModelType::viscType_Newtonian: + get_visc_stress_newt(lDmn.solid_visc.mu, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v); + break; + + case consts::SolidViscosityModelType::viscType_Potential: + get_visc_stress_pot(lDmn.solid_visc.mu, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v); + break; + } +} + }; #endif diff --git a/Code/Source/svFSI/nn.cpp b/Code/Source/svFSI/nn.cpp index d0e5e862..13596cd1 100644 --- a/Code/Source/svFSI/nn.cpp +++ b/Code/Source/svFSI/nn.cpp @@ -439,6 +439,19 @@ void get_xi(const int nsd, consts::ElementType eType, const int eNoN, const Arra xi = xiK; } +/** + * @brief Get shape function gradient w.r.t reference configuration coordinates. + * dN/dX = dN/dxi * dxi/dX, where xi is the parametric coordinates of the parent element + * + * @param eNoN Number of nodes in the element + * @param nsd Number of spatial dimensions + * @param insd + * @param Nxi Shape function gradient w.r.t. parametric coordinates of the parent element + * @param x Element node positions in reference configuration + * @param Nx Shape function gradient w.r.t reference configuration coordinates + * @param Jac Jacobian of element in reference configuration + * @param ks Matrix where each component is (dxi/dX)^2 + */ void gnn(const int eNoN, const int nsd, const int insd, Array& Nxi, Array& x, Array& Nx, double& Jac, Array& ks) { diff --git a/Code/Source/svFSI/read_files.cpp b/Code/Source/svFSI/read_files.cpp index f25d6040..52ae6809 100644 --- a/Code/Source/svFSI/read_files.cpp +++ b/Code/Source/svFSI/read_files.cpp @@ -1299,10 +1299,6 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& } break; - case PhysicalProperyType::solid_viscosity: - rtmp = domain_params->solid_viscosity.value(); - break; - case PhysicalProperyType::source_term: rtmp = domain_params->source_term.value(); break; @@ -1333,7 +1329,13 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType& if ((lEq.dmn[iDmn].phys == EquationType::phys_fluid) || (lEq.dmn[iDmn].phys == EquationType::phys_stokes) || (lEq.dmn[iDmn].phys == EquationType::phys_CMM && !com_mod.cmmInit)) { - read_visc_model(simulation, eq_params, domain_params, lEq.dmn[iDmn]); + read_fluid_visc_model(simulation, eq_params, domain_params, lEq.dmn[iDmn]); + } + + // Set parameters for a fluid viscosity model. + if ((lEq.dmn[iDmn].phys == EquationType::phys_struct) || + (lEq.dmn[iDmn].phys == EquationType::phys_ustruct)) { + read_solid_visc_model(simulation, eq_params, domain_params, lEq.dmn[iDmn]); } } } @@ -2840,11 +2842,11 @@ void read_rmsh(Simulation* simulation, EquationParameters* eq_param) } //----------------- -// read_visc_model +// read_fluid_visc_model //----------------- -// Set the viscosity material model parameters for the given domain. +// Set the fluid viscosity material model parameters for the given domain. // -void read_visc_model(Simulation* simulation, EquationParameters* eq_params, DomainParameters* domain_params, dmnType& lDmn) +void read_fluid_visc_model(Simulation* simulation, EquationParameters* eq_params, DomainParameters* domain_params, dmnType& lDmn) { using namespace consts; @@ -2853,14 +2855,14 @@ void read_visc_model(Simulation* simulation, EquationParameters* eq_params, Doma FluidViscosityModelType vmodel_type; std::string vmodel_str; - if (domain_params->viscosity.model.defined()) { - vmodel_str = domain_params->viscosity.model.value(); + if (domain_params->fluid_viscosity.model.defined()) { + vmodel_str = domain_params->fluid_viscosity.model.value(); std::transform(vmodel_str.begin(), vmodel_str.end(), vmodel_str.begin(), ::tolower); try { vmodel_type = fluid_viscosity_model_name_to_type.at(vmodel_str); } catch (const std::out_of_range& exception) { - throw std::runtime_error("Unknown viscosity model '" + vmodel_str + "'."); + throw std::runtime_error("Unknown fluid viscosity model '" + vmodel_str + "'."); } } else { vmodel_type = FluidViscosityModelType::viscType_Const; @@ -2868,19 +2870,57 @@ void read_visc_model(Simulation* simulation, EquationParameters* eq_params, Doma // Set the parameters for the given viscosity model. // - auto& viscosity_params = domain_params->viscosity; + auto& viscosity_params = domain_params->fluid_viscosity; try { - set_viscosity_props[vmodel_type](simulation, viscosity_params, lDmn); + set_fluid_viscosity_props[vmodel_type](simulation, viscosity_params, lDmn); } catch (const std::bad_function_call& exception) { - throw std::runtime_error("[read_visc_model] Viscosity model '" + vmodel_str + "' is not supported."); + throw std::runtime_error("[read_fluid_visc_model] Viscosity model '" + vmodel_str + "' is not supported."); } - if ((lDmn.phys == EquationType::phys_stokes) && (lDmn.visc.viscType != FluidViscosityModelType::viscType_Const)) { + if ((lDmn.phys == EquationType::phys_stokes) && (lDmn.fluid_visc.viscType != FluidViscosityModelType::viscType_Const)) { throw std::runtime_error("Only constant viscosity is allowed for Stokes flow."); } } +//----------------- +// read_solid_visc_model +//----------------- +// Set the solid viscosity material model parameters for the given domain. +// +void read_solid_visc_model(Simulation* simulation, EquationParameters* eq_params, DomainParameters* domain_params, dmnType& lDmn) +{ + using namespace consts; + + // Get viscosity model. + // + SolidViscosityModelType vmodel_type; + std::string vmodel_str; + + if (domain_params->solid_viscosity.model.defined()) { + vmodel_str = domain_params->solid_viscosity.model.value(); + std::transform(vmodel_str.begin(), vmodel_str.end(), vmodel_str.begin(), ::tolower); + + try { + vmodel_type = solid_viscosity_model_name_to_type.at(vmodel_str); + } catch (const std::out_of_range& exception) { + throw std::runtime_error("Unknown solid viscosity model '" + vmodel_str + "'."); + } + } else { + vmodel_type = SolidViscosityModelType::viscType_Newtonian; + } + + // Set the parameters for the given viscosity model. + // + auto& viscosity_params = domain_params->solid_viscosity; + + try { + set_solid_viscosity_props[vmodel_type](simulation, viscosity_params, lDmn); + } catch (const std::bad_function_call& exception) { + throw std::runtime_error("[read_solid_visc_model] Viscosity model '" + vmodel_str + "' is not supported."); + } +} + //-------------------- // read_wall_props_ff //-------------------- diff --git a/Code/Source/svFSI/read_files.h b/Code/Source/svFSI/read_files.h index 16cdc07a..656b190c 100644 --- a/Code/Source/svFSI/read_files.h +++ b/Code/Source/svFSI/read_files.h @@ -86,7 +86,9 @@ namespace read_files_ns { void read_trac_bcff(ComMod& com_mod, MBType& lMB, faceType& lFa, const std::string& file_name); - void read_visc_model(Simulation* simulation, EquationParameters* eq_params, DomainParameters* domain_params, dmnType& lDmn); + void read_fluid_visc_model(Simulation* simulation, EquationParameters* eq_params, DomainParameters* domain_params, dmnType& lDmn); + + void read_solid_visc_model(Simulation* simulation, EquationParameters* eq_params, DomainParameters* domain_params, dmnType& lDmn); void read_wall_props_ff(ComMod& com_mod, const std::string& file_path, const int iM, const int iFa); diff --git a/Code/Source/svFSI/set_equation_props.h b/Code/Source/svFSI/set_equation_props.h index c9bc5279..1fb5214d 100644 --- a/Code/Source/svFSI/set_equation_props.h +++ b/Code/Source/svFSI/set_equation_props.h @@ -330,11 +330,10 @@ SetEquationPropertiesMapType set_equation_props = { propL[1][n] = PhysicalProperyType::elasticity_modulus; propL[2][n] = PhysicalProperyType::poisson_ratio; propL[3][n] = PhysicalProperyType::damping; - propL[4][n] = PhysicalProperyType::solid_viscosity; - propL[5][n] = PhysicalProperyType::f_x; - propL[6][n] = PhysicalProperyType::f_y; + propL[4][n] = PhysicalProperyType::f_x; + propL[5][n] = PhysicalProperyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[7][n] = PhysicalProperyType::f_z; + propL[6][n] = PhysicalProperyType::f_z; } // Set ustruct properties. @@ -342,13 +341,12 @@ SetEquationPropertiesMapType set_equation_props = { propL[0][n] = PhysicalProperyType::solid_density; propL[1][n] = PhysicalProperyType::elasticity_modulus; propL[2][n] = PhysicalProperyType::poisson_ratio; - propL[3][n] = PhysicalProperyType::solid_viscosity; - propL[4][n] = PhysicalProperyType::ctau_M; - propL[5][n] = PhysicalProperyType::ctau_C; - propL[6][n] = PhysicalProperyType::f_x; - propL[7][n] = PhysicalProperyType::f_y; + propL[3][n] = PhysicalProperyType::ctau_M; + propL[4][n] = PhysicalProperyType::ctau_C; + propL[5][n] = PhysicalProperyType::f_x; + propL[6][n] = PhysicalProperyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[8][n] = PhysicalProperyType::f_z; + propL[7][n] = PhysicalProperyType::f_z; } // Set lElas properties. @@ -570,11 +568,10 @@ SetEquationPropertiesMapType set_equation_props = { propL[1][0] = PhysicalProperyType::damping; propL[2][0] = PhysicalProperyType::elasticity_modulus; propL[3][0] = PhysicalProperyType::poisson_ratio; - propL[4][0] = PhysicalProperyType::solid_viscosity; - propL[5][0] = PhysicalProperyType::f_x; - propL[6][0] = PhysicalProperyType::f_y; + propL[4][0] = PhysicalProperyType::f_x; + propL[5][0] = PhysicalProperyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[7][0] = PhysicalProperyType::f_z; + propL[6][0] = PhysicalProperyType::f_z; } read_domain(simulation, eq_params, lEq, propL); @@ -614,13 +611,12 @@ SetEquationPropertiesMapType set_equation_props = { propL[0][0] = PhysicalProperyType::solid_density; propL[1][0] = PhysicalProperyType::elasticity_modulus; propL[2][0] = PhysicalProperyType::poisson_ratio; - propL[3][0] = PhysicalProperyType::solid_viscosity; - propL[4][0] = PhysicalProperyType::ctau_M; - propL[5][0] = PhysicalProperyType::ctau_C; - propL[6][0] = PhysicalProperyType::f_x; - propL[7][0] = PhysicalProperyType::f_y; + propL[3][0] = PhysicalProperyType::ctau_M; + propL[4][0] = PhysicalProperyType::ctau_C; + propL[5][0] = PhysicalProperyType::f_x; + propL[6][0] = PhysicalProperyType::f_y; if (simulation->com_mod.nsd == 3) { - propL[8][0] = PhysicalProperyType::f_z; + propL[7][0] = PhysicalProperyType::f_z; } read_domain(simulation, eq_params, lEq, propL); diff --git a/Code/Source/svFSI/set_viscosity_props.h b/Code/Source/svFSI/set_viscosity_props.h index 7e03fe68..0730471e 100644 --- a/Code/Source/svFSI/set_viscosity_props.h +++ b/Code/Source/svFSI/set_viscosity_props.h @@ -32,49 +32,51 @@ // S e t V i s c o s i t y P r o p e r t i e s // /////////////////////////////////////////////////////////// // -// The 'set_viscosity_props ' map defined here sets fluid +// The 'set_fluid_viscosity_props ' map defined here sets fluid +// viscosity properties from values read in from a file. +// The 'set_solid_viscosity_props ' map defined here sets solid // viscosity properties from values read in from a file. -using SetViscosityPropertiesMapType = std::map>; -SetViscosityPropertiesMapType set_viscosity_props = { +SetFluidViscosityPropertiesMapType set_fluid_viscosity_props = { //---------------------------// // viscType_Const // //---------------------------// // -{consts::FluidViscosityModelType::viscType_Const, [](Simulation* simulation, ViscosityParameters& params, dmnType& lDmn) -> void +{consts::FluidViscosityModelType::viscType_Const, [](Simulation* simulation, FluidViscosityParameters& params, dmnType& lDmn) -> void { using namespace consts; auto& com_mod = simulation->get_com_mod(); - lDmn.visc.viscType = FluidViscosityModelType::viscType_Const; - lDmn.visc.mu_i = params.newtonian_model.constant_value.value(); + lDmn.fluid_visc.viscType = FluidViscosityModelType::viscType_Const; + lDmn.fluid_visc.mu_i = params.newtonian_model.constant_value.value(); } }, //---------------------------// // viscType_CY // //---------------------------// // -{consts::FluidViscosityModelType::viscType_CY, [](Simulation* simulation, ViscosityParameters& params, dmnType& lDmn) -> void +{consts::FluidViscosityModelType::viscType_CY, [](Simulation* simulation, FluidViscosityParameters& params, dmnType& lDmn) -> void { using namespace consts; auto& com_mod = simulation->get_com_mod(); auto& model_params = params.carreau_yasuda_model; - lDmn.visc.viscType = FluidViscosityModelType::viscType_CY; + lDmn.fluid_visc.viscType = FluidViscosityModelType::viscType_CY; - lDmn.visc.mu_i = model_params.limiting_high_shear_rate_viscosity.value(); - lDmn.visc.mu_o = model_params.limiting_low_shear_rate_viscosity.value(); - lDmn.visc.lam = model_params.shear_rate_tensor_multipler.value(); - lDmn.visc.a = model_params.shear_rate_tensor_exponent.value(); - lDmn.visc.n = model_params.power_law_index.value(); + lDmn.fluid_visc.mu_i = model_params.limiting_high_shear_rate_viscosity.value(); + lDmn.fluid_visc.mu_o = model_params.limiting_low_shear_rate_viscosity.value(); + lDmn.fluid_visc.lam = model_params.shear_rate_tensor_multipler.value(); + lDmn.fluid_visc.a = model_params.shear_rate_tensor_exponent.value(); + lDmn.fluid_visc.n = model_params.power_law_index.value(); - if (lDmn.visc.mu_i > lDmn.visc.mu_o) { + if (lDmn.fluid_visc.mu_i > lDmn.fluid_visc.mu_o) { throw std::runtime_error("Unexpected inputs for Carreau-Yasuda model. " - "High shear-rate viscosity value (" + std::to_string(lDmn.visc.mu_i) + - " should be larger than low shear-rate value " + std::to_string(lDmn.visc.mu_i) + "."); + "High shear-rate viscosity value (" + std::to_string(lDmn.fluid_visc.mu_i) + + " should be larger than low shear-rate value " + std::to_string(lDmn.fluid_visc.mu_i) + "."); } } }, @@ -83,22 +85,61 @@ SetViscosityPropertiesMapType set_viscosity_props = { // viscType_Cass // //---------------------------// // -{consts::FluidViscosityModelType::viscType_Cass, [](Simulation* simulation, ViscosityParameters& params, dmnType& lDmn) -> void +{consts::FluidViscosityModelType::viscType_Cass, [](Simulation* simulation, FluidViscosityParameters& params, dmnType& lDmn) -> void { using namespace consts; auto& com_mod = simulation->get_com_mod(); auto& model_params = params.cassons_model; - lDmn.visc.viscType = FluidViscosityModelType::viscType_Cass; - lDmn.visc.mu_i = model_params.asymptotic_viscosity(); - lDmn.visc.mu_o = model_params.yield_stress(); + lDmn.fluid_visc.viscType = FluidViscosityModelType::viscType_Cass; + lDmn.fluid_visc.mu_i = model_params.asymptotic_viscosity(); + lDmn.fluid_visc.mu_o = model_params.yield_stress(); if (model_params.low_shear_rate_threshold.defined()) { - lDmn.visc.lam = model_params.low_shear_rate_threshold(); + lDmn.fluid_visc.lam = model_params.low_shear_rate_threshold(); } else { - lDmn.visc.lam = 0.5; + lDmn.fluid_visc.lam = 0.5; } } }, }; + + + + + +using SetSolidViscosityPropertiesMapType = std::map>; + +SetSolidViscosityPropertiesMapType set_solid_viscosity_props = { + +//---------------------------// +// viscType_Newtonian // +//---------------------------// +// +{consts::SolidViscosityModelType::viscType_Newtonian, [](Simulation* simulation, SolidViscosityParameters& params, dmnType& lDmn) -> void +{ + using namespace consts; + auto& com_mod = simulation->get_com_mod(); + + lDmn.solid_visc.viscType = SolidViscosityModelType::viscType_Newtonian; + lDmn.solid_visc.mu = params.newtonian_model.constant_value.value(); +} }, + +//---------------------------// +// viscType_Potential // +//---------------------------// +// +{consts::SolidViscosityModelType::viscType_Potential, [](Simulation* simulation, SolidViscosityParameters& params, dmnType& lDmn) -> void +{ + using namespace consts; + auto& com_mod = simulation->get_com_mod(); + + lDmn.solid_visc.viscType = SolidViscosityModelType::viscType_Potential; + lDmn.solid_visc.mu = params.potential_model.constant_value.value(); +} }, + +}; + + diff --git a/Code/Source/svFSI/stokes.cpp b/Code/Source/svFSI/stokes.cpp index 5b9f96ba..a90949cf 100644 --- a/Code/Source/svFSI/stokes.cpp +++ b/Code/Source/svFSI/stokes.cpp @@ -255,7 +255,7 @@ void stokes_2d_c(ComMod& com_mod, const int lStab, const int eNoNw, const int eN int cDmn = com_mod.cDmn; auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; - double mu = dmn.visc.mu_i; + double mu = dmn.fluid_visc.mu_i; double ctM = dmn.prop[PhysicalProperyType::ctau_M]; Vector fb(2); @@ -378,7 +378,7 @@ void stokes_2d_m(ComMod& com_mod, const int eNoNw, const int eNoNq, const double auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; - double mu = dmn.visc.mu_i; + double mu = dmn.fluid_visc.mu_i; Vector fb(2); fb[0] = dmn.prop[PhysicalProperyType::f_x]; fb[1] = dmn.prop[PhysicalProperyType::f_y]; @@ -488,7 +488,7 @@ void stokes_3d_c(ComMod& com_mod, const int lStab, const int eNoNw, const int eN int cDmn = com_mod.cDmn; auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; - double mu = dmn.visc.mu_i; + double mu = dmn.fluid_visc.mu_i; double ctM = dmn.prop[PhysicalProperyType::ctau_M]; Vector fb(3); @@ -604,7 +604,7 @@ void stokes_3d_m(ComMod& com_mod, const int eNoNw, const int eNoNq, const double auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; - double mu = dmn.visc.mu_i; + double mu = dmn.fluid_visc.mu_i; Vector fb(3); fb[0] = dmn.prop[PhysicalProperyType::f_x]; fb[1] = dmn.prop[PhysicalProperyType::f_y]; diff --git a/Code/Source/svFSI/sv_struct.cpp b/Code/Source/svFSI/sv_struct.cpp index f9755d60..224861f1 100644 --- a/Code/Source/svFSI/sv_struct.cpp +++ b/Code/Source/svFSI/sv_struct.cpp @@ -387,7 +387,6 @@ void struct_2d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, // Set parameters // double rho = dmn.prop.at(PhysicalProperyType::solid_density); - double mu = dmn.prop.at(PhysicalProperyType::solid_viscosity); double dmp = dmn.prop.at(PhysicalProperyType::damping); Vector fb({dmn.prop.at(PhysicalProperyType::f_x), dmn.prop.at(PhysicalProperyType::f_y)}); double afu = eq.af * eq.beta*dt*dt; @@ -445,23 +444,17 @@ void struct_2d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, S0(1,0) = S0(0,1); - double Jac = mat_det(F, 2); - auto Fi = mat_inv(F, 2); - - // Viscous contribution - // Velocity gradient in current configuration - auto VxFi = mat_mul(vx, Fi); - - // Deviatoric strain tensor - auto ddev = mat_dev(mat_symm(VxFi,2), 2); - - // 2nd Piola-Kirchhoff stress due to viscosity - auto Svis = mat_mul(ddev, transpose(Fi)); - Svis = 2.0 * mu * Jac * mat_mul(Fi, Svis); - + // 2nd Piola-Kirchhoff stress (S) and material stiffness tensor in Voight notation (Dm) Array S(2,2), Dm(3,3); mat_models::get_pk2cc(com_mod, cep_mod, dmn, F, nFn, fN, ya_g, S, Dm); + // Viscous 2nd Piola-Kirchhoff stress and tangent contributions + Array Svis(2,2); + Array3 Kvis_u(4, eNoN, eNoN); + Array3 Kvis_v(4, eNoN, eNoN); + + mat_models_carray::get_visc_stress_and_tangent<2>(dmn, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v); + // Elastic + Viscous stresses S = S + Svis; @@ -504,22 +497,9 @@ void struct_2d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, Array NxFi(2,eNoN), DdNx(2,eNoN), VxNx(2,eNoN); - for (int a = 0; a < eNoN; a++) { - NxFi(0,a) = Nx(0,a)*Fi(0,0) + Nx(1,a)*Fi(1,0); - NxFi(1,a) = Nx(0,a)*Fi(0,1) + Nx(1,a)*Fi(1,1); - - DdNx(0,a) = ddev(0,0)*NxFi(0,a) + ddev(0,1)*NxFi(1,a); - DdNx(1,a) = ddev(1,0)*NxFi(0,a) + ddev(1,1)*NxFi(1,a); - - VxNx(0,a) = VxFi(0,0)*NxFi(0,a) + VxFi(1,0)*NxFi(1,a); - VxNx(1,a) = VxFi(0,1)*NxFi(0,a) + VxFi(1,1)*NxFi(1,a); - } // Local stiffness tensor - // - double rmu = afu*mu*Jac; - double rmv = afv*mu*Jac; - double T1, Tv, NxNx, NxSNx, BmDBm; + double T1, NxNx, NxSNx, BmDBm; for (int b = 0; b < eNoN; b++) { for (int a = 0; a < eNoN; a++) { @@ -539,51 +519,30 @@ void struct_2d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, DBm(2,0) = Dm(2,0)*Bm(0,0,b) + Dm(2,1)*Bm(1,0,b) + Dm(2,2)*Bm(2,0,b); DBm(2,1) = Dm(2,0)*Bm(0,1,b) + Dm(2,1)*Bm(1,1,b) + Dm(2,2)*Bm(2,1,b); - NxNx = NxFi(0,a)*NxFi(0,b) + NxFi(1,a)*NxFi(1,b); // dM1/du1 // Material stiffness: Bt*D*B BmDBm = Bm(0,0,a)*DBm(0,0) + Bm(1,0,a)*DBm(1,0) + Bm(2,0,a)*DBm(2,0); - // Viscous terms contribution - Tv = (2.0*(DdNx(0,a)*NxFi(0,b) - DdNx(0,b)*NxFi(0,a)) - (NxNx*VxFi(0,0) + NxFi(0,b)*VxNx(0,a) - - NxFi(0,a)*VxNx(0,b))) * rmu + (NxNx) * rmv; - - lK(0,a,b) = lK(0,a,b) + w*(T1 + afu*BmDBm + Tv); + lK(0,a,b) = lK(0,a,b) + w*( T1 + afu*(BmDBm + Kvis_u(0,a,b)) + afv*Kvis_v(0,a,b) ); // dM1/du2 // Material stiffness: Bt*D*B BmDBm = Bm(0,0,a)*DBm(0,1) + Bm(1,0,a)*DBm(1,1) + Bm(2,0,a)*DBm(2,1); - // Viscous terms contribution - Tv = (2.0*(DdNx(0,a)*NxFi(1,b) - DdNx(0,b)*NxFi(1,a)) - - (NxNx*VxFi(0,1) + NxFi(0,b)*VxNx(1,a) - - NxFi(0,a)*VxNx(1,b))) * rmu + (NxFi(1,a)*NxFi(0,b) - - NxFi(0,a)*NxFi(1,b)) * rmv; - - lK(1,a,b) = lK(1,a,b) + w*(afu*BmDBm + Tv); + lK(1,a,b) = lK(1,a,b) + w*( afu*(BmDBm + Kvis_u(1,a,b)) + afv*Kvis_v(1,a,b) ); // dM2/du1 // Material stiffness: Bt*D*B BmDBm = Bm(0,1,a)*DBm(0,0) + Bm(1,1,a)*DBm(1,0) + Bm(2,1,a)*DBm(2,0); - // Viscous terms contribution - Tv = (2.0*(DdNx(1,a)*NxFi(0,b) - DdNx(1,b)*NxFi(0,a)) - - (NxNx*VxFi(1,0) + NxFi(1,b)*VxNx(0,a) - - NxFi(1,a)*VxNx(0,b))) * rmu + (NxFi(0,a)*NxFi(1,b) - - NxFi(1,a)*NxFi(0,b)) * rmv; - - lK(dof+0,a,b) = lK(dof+0,a,b) + w*(afu*BmDBm + Tv); + lK(dof+0,a,b) = lK(dof+0,a,b) + w*( afu*(BmDBm + Kvis_u(dof+0,a,b)) + afv*Kvis_v(dof+0,a,b) ); // dM2/du2 // Material stiffness: Bt*D*B BmDBm = Bm(0,1,a)*DBm(0,1) + Bm(1,1,a)*DBm(1,1) + Bm(2,1,a)*DBm(2,1); - // Viscous terms contribution - Tv = (2.0*(DdNx(1,a)*NxFi(1,b) - DdNx(1,b)*NxFi(1,a)) - (NxNx*VxFi(1,1) + NxFi(1,b)*VxNx(1,a) - - NxFi(1,a)*VxNx(1,b))) * rmu + (NxNx) * rmv; - - lK(dof+1,a,b) = lK(dof+1,a,b) + w*(T1 + afu*BmDBm + Tv); + lK(dof+1,a,b) = lK(dof+1,a,b) + w*( T1 + afu*(BmDBm + Kvis_u(dof+1,a,b)) + afv*Kvis_v(dof+1,a,b) ); } } } @@ -616,7 +575,6 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in // Set parameters // double rho = dmn.prop.at(PhysicalProperyType::solid_density); - double mu = dmn.prop.at(PhysicalProperyType::solid_viscosity); double dmp = dmn.prop.at(PhysicalProperyType::damping); double fb[3]{dmn.prop.at(PhysicalProperyType::f_x), dmn.prop.at(PhysicalProperyType::f_y), @@ -628,7 +586,6 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in #ifdef debug_struct_3d dmsg << "rho: " << rho; - dmsg << "mu: " << mu; dmsg << "dmp: " << dmp; dmsg << "afu: " << afu; dmsg << "afv: " << afv; @@ -656,6 +613,7 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in ud[1] += N(a)*(rho*(al(j,a)-bfl(1,a)) + dmp*yl(j,a)); ud[2] += N(a)*(rho*(al(k,a)-bfl(2,a)) + dmp*yl(k,a)); + // Velocity gradient tensor w.r.t. reference coordinates (dv/dX) vx[0][0] += Nx(0,a)*yl(i,a); vx[0][1] += Nx(1,a)*yl(i,a); vx[0][2] += Nx(2,a)*yl(i,a); @@ -666,6 +624,7 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in vx[2][1] += Nx(1,a)*yl(k,a); vx[2][2] += Nx(2,a)*yl(k,a); + // Deformation gradient tensor (I + du/dX) F[0][0] += Nx(0,a)*dl(i,a); F[0][1] += Nx(1,a)*dl(i,a); F[0][2] += Nx(2,a)*dl(i,a); @@ -691,39 +650,6 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in S0[2][1] = S0[1][2]; S0[0][2] = S0[2][0]; - double Jac = mat_fun_carray::mat_det<3>(F); - - double Fi[3][3]; - mat_fun_carray::mat_inv<3>(F, Fi); - - // Viscous contribution - // Velocity gradient in current configuration - double VxFi[3][3]; - mat_fun_carray::mat_mul(vx, Fi, VxFi); - - // Deviatoric strain tensor - double VxFi_sym[3][3]; - mat_fun_carray::mat_symm<3>(VxFi,VxFi_sym); - - double ddev[3][3]; - mat_fun_carray::mat_dev<3>(VxFi_sym, ddev); - - // 2nd Piola-Kirchhoff stress due to viscosity - double Fi_transp[3][3]; - mat_fun_carray::transpose<3>(Fi, Fi_transp); - - double Svis[3][3]; - mat_fun_carray::mat_mul<3>(ddev, Fi_transp, Svis); - - double Fi_Svis_m[3][3]; - mat_fun_carray::mat_mul<3>(Fi, Svis, Fi_Svis_m); - - for (int i = 0; i < 3; i++) { - for (int j = 0; j < 3; j++) { - Svis[i][j] = 2.0 * mu * Jac * Fi_Svis_m[i][j]; - } - } - // Initialize tensor indexing. mat_fun_carray::ten_init(3); @@ -735,10 +661,18 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in mat_models_carray::get_pk2cc<3>(com_mod, cep_mod, dmn, F, nFn, fN, ya_g, S, Dm); + + // Viscous 2nd Piola-Kirchhoff stress and tangent contributions + Array Svis(3,3); + Array3 Kvis_u(9, eNoN, eNoN); + Array3 Kvis_v(9, eNoN, eNoN); + + mat_models_carray::get_visc_stress_and_tangent<3>(dmn, eNoN, Nx, vx, F, Svis, Kvis_u, Kvis_v); + // Elastic + Viscous stresses for (int i = 0; i < 3; i++) { for (int j = 0; j < 3; j++) { - S[i][j] += Svis[i][j]; + S[i][j] += Svis(i,j); } } @@ -807,31 +741,8 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in Bm(5,2,a) = (Nx(2,a)*F[2][0] + F[2][2]*Nx(0,a)); } - // Below quantities are used for viscous stress contribution - // Shape function gradients in the current configuration - // - Array NxFi(3,eNoN), DdNx(3,eNoN), VxNx(3,eNoN); - - for (int a = 0; a < eNoN; a++) { - NxFi(0,a) = Nx(0,a)*Fi[0][0] + Nx(1,a)*Fi[1][0] + Nx(2,a)*Fi[2][0]; - NxFi(1,a) = Nx(0,a)*Fi[0][1] + Nx(1,a)*Fi[1][1] + Nx(2,a)*Fi[2][1]; - NxFi(2,a) = Nx(0,a)*Fi[0][2] + Nx(1,a)*Fi[1][2] + Nx(2,a)*Fi[2][2]; - - DdNx(0,a) = ddev[0][0]*NxFi(0,a) + ddev[0][1]*NxFi(1,a) + ddev[0][2]*NxFi(2,a); - DdNx(1,a) = ddev[1][0]*NxFi(0,a) + ddev[1][1]*NxFi(1,a) + ddev[1][2]*NxFi(2,a); - DdNx(2,a) = ddev[2][0]*NxFi(0,a) + ddev[2][1]*NxFi(1,a) + ddev[2][2]*NxFi(2,a); - - VxNx(0,a) = VxFi[0][0]*NxFi(0,a) + VxFi[1][0]*NxFi(1,a) + VxFi[2][0]*NxFi(2,a); - VxNx(1,a) = VxFi[0][1]*NxFi(0,a) + VxFi[1][1]*NxFi(1,a) + VxFi[2][1]*NxFi(2,a); - VxNx(2,a) = VxFi[0][2]*NxFi(0,a) + VxFi[1][2]*NxFi(1,a) + VxFi[2][2]*NxFi(2,a); - } - // Local stiffness tensor - double r13 = 1.0 / 3.0; - double r23 = 2.0 / 3.0; - double rmu = afu * mu * Jac; - double rmv = afv * mu * Jac; - double NxSNx, T1, NxNx, BmDBm, Tv; + double NxSNx, T1, NxNx, BmDBm; Array DBm(6,3); @@ -850,7 +761,6 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in // Material Stiffness (Bt*D*B) mat_fun_carray::mat_mul6x3<3>(Dm, Bm.rslice(b), DBm); - NxNx = NxFi(0,a)*NxFi(0,b) + NxFi(1,a)*NxFi(1,b) + NxFi(2,a)*NxFi(2,b); // dM1/du1 // Material stiffness: Bt*D*B @@ -858,11 +768,7 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in Bm(2,0,a)*DBm(2,0) + Bm(3,0,a)*DBm(3,0) + Bm(4,0,a)*DBm(4,0) + Bm(5,0,a)*DBm(5,0); - // Viscous terms contribution - Tv = (2.0*(DdNx(0,a)*NxFi(0,b) - DdNx(0,b)*NxFi(0,a)) - (NxNx*VxFi[0][0] + NxFi(0,b)*VxNx(0,a) - - r23*NxFi(0,a)*VxNx(0,b))) * rmu + (r13*NxFi(0,a)*NxFi(0,b) + NxNx) * rmv; - - lK(0,a,b) = lK(0,a,b) + w*(T1 + afu*BmDBm + Tv); + lK(0,a,b) = lK(0,a,b) + w*( T1 + afu*(BmDBm + Kvis_u(0,a,b)) + afv*Kvis_v(0,a,b) ); // dM1/du2 // Material stiffness: Bt*D*B @@ -870,13 +776,8 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in Bm(2,0,a)*DBm(2,1) + Bm(3,0,a)*DBm(3,1) + Bm(4,0,a)*DBm(4,1) + Bm(5,0,a)*DBm(5,1); - // Viscous terms contribution - Tv = (2.0*(DdNx(0,a)*NxFi(1,b) - DdNx(0,b)*NxFi(1,a)) - - (NxNx*VxFi[0][1] + NxFi(0,b)*VxNx(1,a) - - r23*NxFi(0,a)*VxNx(1,b))) * rmu - + (NxFi(1,a)*NxFi(0,b) - r23*NxFi(0,a)*NxFi(1,b)) * rmv; - lK(1,a,b) = lK(1,a,b) + w*(afu*BmDBm + Tv); + lK(1,a,b) = lK(1,a,b) + w*( afu*(BmDBm + Kvis_u(1,a,b)) + afv*(Kvis_v(1,a,b)) ); // dM1/du3 // Material stiffness: Bt*D*B @@ -884,13 +785,8 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in Bm(2,0,a)*DBm(2,2) + Bm(3,0,a)*DBm(3,2) + Bm(4,0,a)*DBm(4,2) + Bm(5,0,a)*DBm(5,2); - // Viscous terms contribution - Tv = (2.0*(DdNx(0,a)*NxFi(2,b) - DdNx(0,b)*NxFi(2,a)) - - (NxNx*VxFi[0][2] + NxFi(0,b)*VxNx(2,a) - - r23*NxFi(0,a)*VxNx(2,b))) * rmu + - (NxFi(2,a)*NxFi(0,b) - r23*NxFi(0,a)*NxFi(2,b)) * rmv; - lK(2,a,b) = lK(2,a,b) + w*(afu*BmDBm + Tv); + lK(2,a,b) = lK(2,a,b) + w*( afu*(BmDBm + Kvis_u(2,a,b)) + afv*Kvis_v(2,a,b) ); // dM2/du1 // Material stiffness: Bt*D*B @@ -898,13 +794,8 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in Bm(2,1,a)*DBm(2,0) + Bm(3,1,a)*DBm(3,0) + Bm(4,1,a)*DBm(4,0) + Bm(5,1,a)*DBm(5,0); - // Viscous terms contribution - Tv = (2.0*(DdNx(1,a)*NxFi(0,b) - DdNx(1,b)*NxFi(0,a)) - - (NxNx*VxFi[1][0] + NxFi(1,b)*VxNx(0,a) - - r23*NxFi(1,a)*VxNx(0,b))) * rmu + - (NxFi(0,a)*NxFi(1,b) - r23*NxFi(1,a)*NxFi(0,b)) * rmv; - lK(dof+0,a,b) = lK(dof+0,a,b) + w*(afu*BmDBm + Tv); + lK(dof+0,a,b) = lK(dof+0,a,b) + w*( afu*(BmDBm + Kvis_u(dof+0,a,b)) + afv*Kvis_v(dof+0,a,b) ); // dM2/du2 // Material stiffness: Bt*D*B @@ -912,13 +803,8 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in Bm(2,1,a)*DBm(2,1) + Bm(3,1,a)*DBm(3,1) + Bm(4,1,a)*DBm(4,1) + Bm(5,1,a)*DBm(5,1); - // Viscous terms contribution - Tv = (2.0*(DdNx(1,a)*NxFi(1,b) - DdNx(1,b)*NxFi(1,a)) - - (NxNx*VxFi[1][1] + NxFi(1,b)*VxNx(1,a) - - r23*NxFi(1,a)*VxNx(1,b))) * rmu + - (r13*NxFi(1,a)*NxFi(1,b) + NxNx) * rmv; - lK(dof+1,a,b) = lK(dof+1,a,b) + w*(T1 + afu*BmDBm + Tv); + lK(dof+1,a,b) = lK(dof+1,a,b) + w*(T1 + afu*(BmDBm + Kvis_u(dof+1,a,b)) + afv*Kvis_v(dof+1,a,b) ); // dM2/du3 // Material stiffness: Bt*D*B @@ -926,13 +812,8 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in Bm(2,1,a)*DBm(2,2) + Bm(3,1,a)*DBm(3,2) + Bm(4,1,a)*DBm(4,2) + Bm(5,1,a)*DBm(5,2); - // Viscous terms contribution - Tv = (2.0*(DdNx(1,a)*NxFi(2,b) - DdNx(1,b)*NxFi(2,a)) - - (NxNx*VxFi[1][2] + NxFi(1,b)*VxNx(2,a) - - r23*NxFi(1,a)*VxNx(2,b))) * rmu + (NxFi(2,a)*NxFi(1,b) - - r23*NxFi(1,a)*NxFi(2,b)) * rmv; - lK(dof+2,a,b) = lK(dof+2,a,b) + w*(afu*BmDBm + Tv); + lK(dof+2,a,b) = lK(dof+2,a,b) + w*( afu*(BmDBm + Kvis_u(dof+2,a,b)) + afv*Kvis_v(dof+2,a,b) ); // dM3/du1 // Material stiffness: Bt*D*B @@ -940,13 +821,8 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in Bm(2,2,a)*DBm(2,0) + Bm(3,2,a)*DBm(3,0) + Bm(4,2,a)*DBm(4,0) + Bm(5,2,a)*DBm(5,0); - // Viscous terms contribution - Tv = (2.0*(DdNx(2,a)*NxFi(0,b) - DdNx(2,b)*NxFi(0,a)) - - (NxNx*VxFi[2][0] + NxFi(2,b)*VxNx(0,a) - - r23*NxFi(2,a)*VxNx(0,b))) * rmu + (NxFi(0,a)*NxFi(2,b) - - r23*NxFi(2,a)*NxFi(0,b)) * rmv; - lK(2*dof+0,a,b) = lK(2*dof+0,a,b) + w*(afu*BmDBm + Tv); + lK(2*dof+0,a,b) = lK(2*dof+0,a,b) + w*( afu*(BmDBm + Kvis_u(2*dof+0,a,b)) + afv*Kvis_v(2*dof+0,a,b) ); // dM3/du2 // Material stiffness: Bt*D*B @@ -954,13 +830,8 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in Bm(2,2,a)*DBm(2,1) + Bm(3,2,a)*DBm(3,1) + Bm(4,2,a)*DBm(4,1) + Bm(5,2,a)*DBm(5,1); - // Viscous terms contribution - Tv = (2.0*(DdNx(2,a)*NxFi(1,b) - DdNx(2,b)*NxFi(1,a)) - - (NxNx*VxFi[2][1] + NxFi(2,b)*VxNx(1,a) - - r23*NxFi(2,a)*VxNx(1,b))) * rmu + (NxFi(1,a)*NxFi(2,b) - - r23*NxFi(2,a)*NxFi(1,b)) * rmv; - lK(2*dof+1,a,b) = lK(2*dof+1,a,b) + w*(afu*BmDBm + Tv); + lK(2*dof+1,a,b) = lK(2*dof+1,a,b) + w*( afu*(BmDBm + Kvis_u(2*dof+1,a,b)) + afv*Kvis_v(2*dof+1,a,b) ); // dM3/du3 // Material stiffness: Bt*D*B @@ -968,13 +839,7 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in Bm(2,2,a)*DBm(2,2) + Bm(3,2,a)*DBm(3,2) + Bm(4,2,a)*DBm(4,2) + Bm(5,2,a)*DBm(5,2); - // Viscous terms contribution - Tv = (2.0*(DdNx(2,a)*NxFi(2,b) - DdNx(2,b)*NxFi(2,a)) - - (NxNx*VxFi[2][2] + NxFi(2,b)*VxNx(2,a) - - r23*NxFi(2,a)*VxNx(2,b))) * rmu + - (r13*NxFi(2,a)*NxFi(2,b) + NxNx) * rmv; - - lK(2*dof+2,a,b) = lK(2*dof+2,a,b) + w*(T1 + afu*BmDBm + Tv); + lK(2*dof+2,a,b) = lK(2*dof+2,a,b) + w*( T1 + afu*(BmDBm + Kvis_u(2*dof+2,a,b)) + afv*Kvis_v(2*dof+2,a,b) ); } } @@ -982,6 +847,7 @@ void struct_3d_carray(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const in /// @brief Reproduces Fortran 'STRUCT3D' subroutine. // +// DEPRECATED: This function is deprecated and will be removed in the future. void struct_3d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, const double w, const Vector& N, const Array& Nx, const Array& al, const Array& yl, const Array& dl, const Array& bfl, const Array& fN, const Array& pS0l, @@ -1009,7 +875,7 @@ void struct_3d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, // Set parameters // double rho = dmn.prop.at(PhysicalProperyType::solid_density); - double mu = dmn.prop.at(PhysicalProperyType::solid_viscosity); + double mu = 0.0; double dmp = dmn.prop.at(PhysicalProperyType::damping); Vector fb({dmn.prop.at(PhysicalProperyType::f_x), dmn.prop.at(PhysicalProperyType::f_y), @@ -1021,7 +887,6 @@ void struct_3d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, #ifdef debug_struct_3d dmsg << "rho: " << rho; - dmsg << "mu: " << mu; dmsg << "dmp: " << dmp; dmsg << "afu: " << afu; dmsg << "afv: " << afv; diff --git a/Code/Source/svFSI/ustruct.cpp b/Code/Source/svFSI/ustruct.cpp index 25ab4162..cf8602f9 100644 --- a/Code/Source/svFSI/ustruct.cpp +++ b/Code/Source/svFSI/ustruct.cpp @@ -52,6 +52,7 @@ #include "fs.h" #include "mat_fun.h" #include "mat_models.h" +#include "mat_models_carray.h" #include "nn.h" #include "utils.h" @@ -384,13 +385,13 @@ void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const auto N0 = fs[0].N.col(g); auto N1 = fs[1].N.col(g); ustruct_3d_c(com_mod, cep_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, Jac, N0, N1, Nwx, - Nqx, al, yl, dl, bfl, ksix, lR, lK, lKd); + Nqx, al, yl, dl, bfl, lR, lK, lKd); } else if (nsd == 2) { auto N0 = fs[0].N.col(g); auto N1 = fs[1].N.col(g); ustruct_2d_c(com_mod, cep_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, Jac, N0, N1, Nwx, - Nqx, al, yl, dl, bfl, ksix, lR, lK, lKd); + Nqx, al, yl, dl, bfl, lR, lK, lKd); } } // for g = 0 to fs[1].nG @@ -427,7 +428,7 @@ int get_col_ptr(ComMod& com_mod, const int rowN, const int colN) void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, const double w, const double Je, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& al, const Array& yl, - const Array& dl, const Array& bfl, const Array& Kxi, Array& lR, Array3& lK, + const Array& dl, const Array& bfl, Array& lR, Array3& lK, Array3& lKd) { using namespace consts; @@ -449,8 +450,6 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in int cDmn = com_mod.cDmn; auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; - double ctV = 36.0; - double mu = dmn.prop[PhysicalProperyType::solid_viscosity]; Vector fb(2); fb[0] = dmn.prop[PhysicalProperyType::f_x]; @@ -527,21 +526,10 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in // double tauM = 0.0; double tauC = 0.0; - Array Kx(2,2); if (vmsFlag) { mat_models::get_tau(com_mod, eq.dmn[cDmn], Jac, Je, tauM, tauC); - // Stabilization parameter for solid viscosity - Kx = mat_mul(Kxi, Fi); - Kx = mat_mul(transpose(Fi), Kx); - - double tauV = Kx(0,0)*Kx(0,0) + Kx(1,0)*Kx(1,0) + - Kx(0,1)*Kx(0,1) + Kx(1,1)*Kx(1,1); - - tauV = ctV * mu * mu * tauV; - tauM = 1.0 / sqrt( tauV + pow(1.0/tauM,2.0)); - } else { tauM = 0.0; tauC = 0.0; @@ -644,7 +632,7 @@ void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, const double w, const double Je, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& al, const Array& yl, - const Array& dl, const Array& bfl, const Array& Kxi, Array& lR, Array3& lK, + const Array& dl, const Array& bfl, Array& lR, Array3& lK, Array3& lKd) { using namespace consts; @@ -666,8 +654,6 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in int cDmn = com_mod.cDmn; auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; - double ctV = 36.0; - double mu = dmn.prop[PhysicalProperyType::solid_viscosity]; Vector fb(3); fb[0] = dmn.prop[PhysicalProperyType::f_x]; @@ -763,22 +749,10 @@ void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in // double tauM = 0.0; double tauC = 0.0; - Array Kx(3,3); if (vmsFlag) { mat_models::get_tau(com_mod, eq.dmn[cDmn], Jac, Je, tauM, tauC); - // Stabilization parameter for solid viscosity - Kx = mat_mul(Kxi, Fi); - Kx = mat_mul(transpose(Fi), Kx); - - double tauV = Kx(0,0)*Kx(0,0) + Kx(1,0)*Kx(1,0) + Kx(2,0)*Kx(2,0) + - Kx(0,1)*Kx(0,1) + Kx(1,1)*Kx(1,1) + Kx(2,1)*Kx(2,1) + - Kx(0,2)*Kx(0,2) + Kx(1,2)*Kx(1,2) + Kx(2,2)*Kx(2,2); - - tauV = ctV * mu * mu * tauV; - tauM = 1.0 / sqrt( tauV + pow(1.0/tauM,2.0)); - } else { tauM = 0.0; tauC = 0.0; @@ -931,7 +905,6 @@ void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in fb[0] = dmn.prop[PhysicalProperyType::f_x]; fb[1] = dmn.prop[PhysicalProperyType::f_y]; - double mu = dmn.prop[PhysicalProperyType::solid_viscosity]; double am = eq.am; double af = eq.af * eq.gam * dt; double afm = af / am; @@ -997,16 +970,12 @@ void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in double Ja = 0; mat_models::get_pk2cc_dev(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, fN, ya_g, Siso, Dm, Ja); - // Viscous contribution - // Velocity gradient in current configuration - auto VxFi = mat_mul(vx, Fi); - - // Deviatoric strain tensor - auto ddev = mat_dev(mat_symm(VxFi,2), 2); - - // 2nd Piola-Kirchhoff stress due to viscosity - auto Svis = mat_mul(ddev, transpose(Fi)); - Svis = 2.0*mu*Jac*mat_mul(Fi, Svis); + // Viscous 2nd Piola-Kirchhoff stress and tangent contributions + Array Svis(2,2); + Array3 Kvis_u(4, eNoNw, eNoNw); + Array3 Kvis_v(4, eNoNw, eNoNw); + + mat_models_carray::get_visc_stress_and_tangent<2>(dmn, eNoNw, Nwx, vx, F, Svis, Kvis_u, Kvis_v); // Compute rho and beta depending on the volumetric penalty model // @@ -1027,7 +996,7 @@ void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in tauC = 0.0; } - // Total isochoric 2nd Piola-Kirchhoff stress + // Total isochoric 2nd Piola-Kirchhoff stress (Elastic + Viscous) Siso = Siso + Svis; // Deviatoric 1st Piola-Kirchhoff tensor (P) @@ -1044,6 +1013,8 @@ void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in NxFi(1,a) = Nwx(0,a)*Fi(0,1) + Nwx(1,a)*Fi(1,1); } + // Velocity gradient in current configuration + auto VxFi = mat_mul(vx, Fi); double rC = beta*pd + VxFi(1,1) + VxFi(2,2); double rCl = -p + tauC*rC; @@ -1076,19 +1047,16 @@ void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in Bm(2,1,a) = Nwx(2,a)*F(1,2) + F(1,0)*Nwx(1,a); } - Array VxNx(2,eNoNw), DdNx(2,eNoNw); + Array VxNx(2,eNoNw); for (int a = 0; a < eNoNw; a++) { - DdNx(0,a) = ddev(0,0)*NxFi(0,a) + ddev(0,1)*NxFi(1,a); - DdNx(1,a) = ddev(1,0)*NxFi(0,a) + ddev(1,1)*NxFi(1,a); - VxNx(0,a) = VxFi(0,0)*NxFi(0,a) + VxFi(1,0)*NxFi(1,a); VxNx(1,a) = VxFi(0,1)*NxFi(0,a) + VxFi(1,1)*NxFi(1,a); } // Tangent (stiffness) matrices // - double NxSNx{0.0}, NxNx{0.0}, BtDB{0.0}, + double NxSNx{0.0}, BtDB{0.0}, T1{0.0}, T2{0.0}, T3{0.0}, Tv{0.0}, Ku{0.0}; @@ -1110,23 +1078,19 @@ void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in DBm(2,0) = Dm(2,0)*Bm(0,0,b) + Dm(2,1)*Bm(1,0,b) + Dm(2,2)*Bm(2,0,b); DBm(2,1) = Dm(2,0)*Bm(0,1,b) + Dm(2,1)*Bm(1,1,b) + Dm(2,2)*Bm(2,1,b); - NxNx = NxFi(0,a)*NxFi(0,b) + NxFi(1,a)*NxFi(1,b); // dM1_dV1 + af/am *dM_1/dU_1 // BtDB = Bm(0,0,a)*DBm(0,0) + Bm(1,0,a)*DBm(1,0) + Bm(2,0,a)*DBm(2,0); T1 = Jac*rho*vd(0)*Nw(a)*NxFi(0,b); T2 = -tauC*Jac*NxFi(0,a)*VxNx(0,b); - Tv = (1.0*(DdNx(0,a)*NxFi(0,b) - DdNx(0,b)*NxFi(0,a)) - - (NxNx*VxFi(0,0) + NxFi(0,b)*VxNx(0,a) - - NxFi(0,a)*VxNx(0,b)))*mu*Jac; - Ku = w*af*(T1 + T2 + Tv + BtDB + NxSNx); + Ku = w*af*(T1 + T2 + BtDB + NxSNx + Kvis_u(0,a,b)); lKd(0,a,b) = lKd(0,a,b) + Ku; T1 = am*Jac*rho*Nw(a)*Nw(b); T2 = T1 + af*Jac*tauC*rho*NxFi(0,a)*NxFi(0,b); - Tv = af*mu*Jac*NxNx; + Tv = af*Kvis_v(0,a,b); lK(0,a,b) = lK(0,a,b) + w*(T2 + Tv) + afm*Ku; // dM_1/dV_2 + af/am *dM_1/dU_2 @@ -1135,16 +1099,12 @@ void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T1 = Jac*rho*vd(0)*Nw(a)*NxFi(1,b); T2 = -tauC*Jac*NxFi(0,a)*VxNx(1,b); T3 = Jac*rCl*(NxFi(0,a)*NxFi(1,b) - NxFi(1,a)*NxFi(0,b)); - Tv = (1.0*(DdNx(0,a)*NxFi(1,b) - DdNx(0,b)*NxFi(1,a)) - - (NxNx*VxFi(0,1) + NxFi(0,b)*VxNx(1,a) - - NxFi(0,a)*VxNx(1,b)))*mu*Jac; - Ku = w*af*(T1 + T2 + T3 + Tv + BtDB); + Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(1,a,b)); lKd(1,a,b) = lKd(1,a,b) + Ku; T2 = af*Jac*tauC*rho*NxFi(0,a)*NxFi(1,b); - Tv = af*mu*Jac*(NxFi(1,a)*NxFi(0,b) - - NxFi(0,a)*NxFi(1,b)); + Tv = af*Kvis_v(1,a,b); lK(1,a,b) = lK(1,a,b) + w*(T2 + Tv) + afm*Ku; // dM_2/dV_1 + af/am *dM_2/dU_1 @@ -1153,16 +1113,12 @@ void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T1 = Jac*rho*vd(1)*Nw(a)*NxFi(0,b); T2 = -tauC*Jac*NxFi(1,a)*VxNx(0,b); T3 = Jac*rCl*(NxFi(1,a)*NxFi(0,b) - NxFi(0,a)*NxFi(1,b)); - Tv = (1.0*(DdNx(1,a)*NxFi(0,b) - DdNx(1,b)*NxFi(0,a)) - - (NxNx*VxFi(1,0) + NxFi(1,b)*VxNx(0,a) - - NxFi(1,a)*VxNx(0,b)))*mu*Jac; - Ku = w*af*(T1 + T2 + T3 + Tv + BtDB); + Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(2,a,b)); lKd(2,a,b) = lKd(2,a,b) + Ku; T2 = af*Jac*tauC*rho*NxFi(1,a)*NxFi(0,b); - Tv = af*mu*Jac*(NxFi(0,a)*NxFi(1,b) - - NxFi(1,a)*NxFi(0,b)); + Tv = af*Kvis_v(2,a,b); lK(3,a,b) = lK(3,a,b) + w*(T2 + Tv) + afm*Ku; // dM_2/dV_2 + af/am *dM_2/dU_2 @@ -1170,16 +1126,13 @@ void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in BtDB = Bm(0,1,a)*DBm(0,1) + Bm(1,1,a)*DBm(1,1) + Bm(2,1,a)*DBm(2,1); T1 = Jac*rho*vd(1)*Nw(a)*NxFi(1,b); T2 = -tauC*Jac*NxFi(1,a)*VxNx(1,b); - Tv = (1.0*(DdNx(1,a)*NxFi(1,b) - DdNx(1,b)*NxFi(1,a)) - - (NxNx*VxFi(1,1) + NxFi(1,b)*VxNx(1,a) - - NxFi(1,a)*VxNx(1,b)))*mu*Jac; - Ku = w*af*(T1 + T2 + Tv + BtDB + NxSNx); + Ku = w*af*(T1 + T2 + BtDB + NxSNx + Kvis_u(3,a,b)); lKd(3,a,b) = lKd(3,a,b) + Ku; T1 = am*Jac*rho*Nw(a)*Nw(b); T2 = T1 + af*Jac*tauC*rho*NxFi(1,a)*NxFi(1,b); - Tv = af*mu*Jac*NxNx; + Tv = af*Kvis_v(3,a,b); lK(4,a,b) = lK(4,a,b) + w*(T2 + Tv) + afm*Ku; } } @@ -1227,8 +1180,6 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in const double dt = com_mod.dt; // Define parameters - // - double mu = dmn.prop[PhysicalProperyType::solid_viscosity]; Vector fb(3); fb[0] = dmn.prop[PhysicalProperyType::f_x]; @@ -1320,17 +1271,13 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in double Ja = 0; mat_models::get_pk2cc_dev(com_mod, cep_mod, eq.dmn[cDmn], F, nFn, fN, ya_g, Siso, Dm, Ja); - // Viscous contribution - // - // Velocity gradient in current configuration - auto VxFi = mat_mul(vx, Fi); - - // Deviatoric strain tensor - auto ddev = mat_dev(mat_symm(VxFi,3), 3); + // Viscous 2nd Piola-Kirchhoff stress and tangent contributions + Array Svis(3,3); + Array3 Kvis_u(9, eNoNw, eNoNw); + Array3 Kvis_v(9, eNoNw, eNoNw); + + mat_models_carray::get_visc_stress_and_tangent<3>(dmn, eNoNw, Nwx, vx, F, Svis, Kvis_u, Kvis_v); - // 2nd Piola-Kirchhoff stress due to viscosity - auto Svis = mat_mul(ddev, transpose(Fi)); - Svis = 2.0 * mu * Jac * mat_mul(Fi, Svis); // Compute rho and beta depending on the volumetric penalty model // @@ -1351,7 +1298,7 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in tauC = 0.0; } - // Total isochoric 2nd Piola-Kirchhoff stress + // Total isochoric 2nd Piola-Kirchhoff stress (Elastic + Viscous) Siso = Siso + Svis; // Deviatoric 1st Piola-Kirchhoff tensor (P) @@ -1368,6 +1315,8 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in NxFi(2,a) = Nwx(0,a)*Fi(0,2) + Nwx(1,a)*Fi(1,2) + Nwx(2,a)*Fi(2,2); } + // Velocity gradient in current configuration + auto VxFi = mat_mul(vx, Fi); double rC = beta*pd + VxFi(0,0) + VxFi(1,1) + VxFi(2,2); double rCl = -p + tauC*rC; @@ -1422,13 +1371,9 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in Bm(5,2,a) = (Nwx(2,a)*F(2,0) + F(2,2)*Nwx(0,a)); } - Array VxNx(3,eNoNw), DdNx(3,eNoNw); + Array VxNx(3,eNoNw); for (int a = 0; a < eNoNw; a++) { - DdNx(0,a) = ddev(0,0)*NxFi(0,a) + ddev(0,1)*NxFi(1,a) + ddev(0,2)*NxFi(2,a); - DdNx(1,a) = ddev(1,0)*NxFi(0,a) + ddev(1,1)*NxFi(1,a) + ddev(1,2)*NxFi(2,a); - DdNx(2,a) = ddev(2,0)*NxFi(0,a) + ddev(2,1)*NxFi(1,a) + ddev(2,2)*NxFi(2,a); - VxNx(0,a) = VxFi(0,0)*NxFi(0,a) + VxFi(1,0)*NxFi(1,a) + VxFi(2,0)*NxFi(2,a); VxNx(1,a) = VxFi(0,1)*NxFi(0,a) + VxFi(1,1)*NxFi(1,a) + VxFi(2,1)*NxFi(2,a); VxNx(2,a) = VxFi(0,2)*NxFi(0,a) + VxFi(1,2)*NxFi(1,a) + VxFi(2,2)*NxFi(2,a); @@ -1438,7 +1383,7 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in // double r13 = 1.0 / 3.0; double r23 = 2.0 / 3.0; - double NxSNx{0.0}, BtDB{0.0}, NxNx{0.0}; + double NxSNx{0.0}, BtDB{0.0}; double Tv{0.0}, Ku{0.0}; for (int b = 0; b < eNoNw; b++) { @@ -1450,7 +1395,6 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in + Nwx(2,a)*Siso(2,1)*Nwx(1,b) + Nwx(2,a)*Siso(2,2)*Nwx(2,b); auto DBm = mat_mul(Dm, Bm.rslice(b)); - NxNx = NxFi(0,a)*NxFi(0,b) + NxFi(1,a)*NxFi(1,b) + NxFi(2,a)*NxFi(2,b); // dM1_dV1 + af/am *dM_1/dU_1 BtDB = Bm(0,0,a)*DBm(0,0) + Bm(1,0,a)*DBm(1,0) + @@ -1458,16 +1402,13 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in Bm(4,0,a)*DBm(4,0) + Bm(5,0,a)*DBm(5,0); T1 = Jac*rho*vd(0)*Nw(a)*NxFi(0,b); T2 = -tauC*Jac*NxFi(0,a)*VxNx(0,b); - Tv = (2.0*(DdNx(0,a)*NxFi(0,b) - DdNx(0,b)*NxFi(0,a)) - - (NxNx*VxFi(0,0) + NxFi(0,b)*VxNx(0,a) - - r23*NxFi(0,a)*VxNx(0,b)))*mu*Jac; - Ku = w*af*(T1 + T2 + Tv + BtDB + NxSNx); + Ku = w*af*(T1 + T2 + BtDB + NxSNx + Kvis_u(0,a,b)); lKd(0,a,b) = lKd(0,a,b) + Ku; T1 = am*Jac*rho*Nw(a)*Nw(b); T2 = T1 + af*Jac*tauC*rho*NxFi(0,a)*NxFi(0,b); - Tv = af*mu*Jac*(r13*NxFi(0,a)*NxFi(0,b) + NxNx); + Tv = af*Kvis_v(0,a,b); lK(0,a,b) = lK(0,a,b) + w*(T2 + Tv) + afm*Ku; // dM_1/dV_2 + af/am *dM_1/dU_2 @@ -1477,16 +1418,12 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T1 = Jac*rho*vd(0)*Nw(a)*NxFi(1,b); T2 = -tauC*Jac*NxFi(0,a)*VxNx(1,b); T3 = Jac*rCl*(NxFi(0,a)*NxFi(1,b) - NxFi(1,a)*NxFi(0,b)); - Tv = (2.0*(DdNx(0,a)*NxFi(1,b) - DdNx(0,b)*NxFi(1,a)) - - (NxNx*VxFi(0,1) + NxFi(0,b)*VxNx(1,a) - - r23*NxFi(0,a)*VxNx(1,b)))*mu*Jac; - Ku = w*af*(T1 + T2 + T3 + Tv + BtDB); + Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(1,a,b)); lKd(1,a,b) = lKd(1,a,b) + Ku; T2 = af*Jac*tauC*rho*NxFi(0,a)*NxFi(1,b); - Tv = af*mu*Jac*(NxFi(1,a)*NxFi(0,b) - - r23*NxFi(0,a)*NxFi(1,b)); + Tv = af*Kvis_v(1,a,b); lK(1,a,b) = lK(1,a,b) + w*(T2 + Tv) + afm*Ku; // dM_1/dV_3 + af/am *dM_1/dU_3 @@ -1497,16 +1434,12 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T1 = Jac*rho*vd(0)*Nw(a)*NxFi(2,b); T2 = -tauC*Jac*NxFi(0,a)*VxNx(2,b); T3 = Jac*rCl*(NxFi(0,a)*NxFi(2,b) - NxFi(2,a)*NxFi(0,b)); - Tv = (2.0*(DdNx(0,a)*NxFi(2,b) - DdNx(0,b)*NxFi(2,a)) - - (NxNx*VxFi(0,2) + NxFi(0,b)*VxNx(2,a) - - r23*NxFi(0,a)*VxNx(2,b)))*mu*Jac; - Ku = w*af*(T1 + T2 + T3 + Tv + BtDB); + Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(2,a,b)); lKd(2,a,b) = lKd(2,a,b) + Ku; T2 = af*Jac*tauC*rho*NxFi(0,a)*NxFi(2,b); - Tv = af*mu*Jac*(NxFi(2,a)*NxFi(0,b) - - r23*NxFi(0,a)*NxFi(2,b)); + Tv = af*Kvis_v(2,a,b); lK(2,a,b) = lK(2,a,b) + w*(T2 + Tv) + afm*Ku; // dM_2/dV_1 + af/am *dM_2/dU_1 @@ -1518,16 +1451,12 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T1 = Jac*rho*vd(1)*Nw(a)*NxFi(0,b); T2 = -tauC*Jac*NxFi(1,a)*VxNx(0,b); T3 = Jac*rCl*(NxFi(1,a)*NxFi(0,b) - NxFi(0,a)*NxFi(1,b)); - Tv = (2.0*(DdNx(1,a)*NxFi(0,b) - DdNx(1,b)*NxFi(0,a)) - - (NxNx*VxFi(1,0) + NxFi(1,b)*VxNx(0,a) - - r23*NxFi(1,a)*VxNx(0,b)))*mu*Jac; - Ku = w*af*(T1 + T2 + T3 + Tv + BtDB); + Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(3,a,b)); lKd(3,a,b) = lKd(3,a,b) + Ku; T2 = af*Jac*tauC*rho*NxFi(1,a)*NxFi(0,b); - Tv = af*mu*Jac*(NxFi(0,a)*NxFi(1,b) - - r23*NxFi(1,a)*NxFi(0,b)); + Tv = af*Kvis_v(3,a,b); lK(4,a,b) = lK(4,a,b) + w*(T2 + Tv) + afm*Ku; @@ -1541,16 +1470,13 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T2 = -tauC*Jac*NxFi(1,a)*VxNx(1,b); - Tv = (2.0*(DdNx(1,a)*NxFi(1,b) - DdNx(1,b)*NxFi(1,a)) - - (NxNx*VxFi(1,1) + NxFi(1,b)*VxNx(1,a) - - r23*NxFi(1,a)*VxNx(1,b)))*mu*Jac; - Ku = w*af*(T1 + T2 + Tv + BtDB + NxSNx); + Ku = w*af*(T1 + T2 + BtDB + NxSNx + Kvis_u(4,a,b)); lKd(4,a,b) = lKd(4,a,b) + Ku; T1 = am*Jac*rho*Nw(a)*Nw(b); T2 = T1 + af*Jac*tauC*rho*NxFi(1,a)*NxFi(1,b); - Tv = af*mu*Jac*(r13*NxFi(1,a)*NxFi(1,b) + NxNx); + Tv = af*Kvis_v(4,a,b); lK(5,a,b) = lK(5,a,b) + w*(T2 + Tv) + afm*Ku; @@ -1564,16 +1490,12 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T2 = -tauC*Jac*NxFi(1,a)*VxNx(2,b); T3 = Jac*rCl*(NxFi(1,a)*NxFi(2,b) - NxFi(2,a)*NxFi(1,b)); - Tv = (2.0*(DdNx(1,a)*NxFi(2,b) - DdNx(1,b)*NxFi(2,a)) - - (NxNx*VxFi(1,2) + NxFi(1,b)*VxNx(2,a) - - r23*NxFi(1,a)*VxNx(2,b)))*mu*Jac; - Ku = w*af*(T1 + T2 + T3 + Tv + BtDB); + Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(5,a,b)); lKd(5,a,b) = lKd(5,a,b) + Ku; T2 = af*Jac*tauC*rho*NxFi(1,a)*NxFi(2,b); - Tv = af*mu*Jac*(NxFi(2,a)*NxFi(1,b) - - r23*NxFi(1,a)*NxFi(2,b)); + Tv = af*Kvis_v(5,a,b); lK(6,a,b) = lK(6,a,b) + w*(T2 + Tv) + afm*Ku; // dM_3/dV_1 + af/am *dM_3/dU_1 @@ -1585,17 +1507,12 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T1 = Jac*rho*vd(2)*Nw(a)*NxFi(0,b); T2 = -tauC*Jac*NxFi(2,a)*VxNx(0,b); T3 = Jac*rCl*(NxFi(2,a)*NxFi(0,b) - NxFi(0,a)*NxFi(2,b)); - - Tv = (2.0*(DdNx(2,a)*NxFi(0,b) - DdNx(2,b)*NxFi(0,a)) - - (NxNx*VxFi(2,0) + NxFi(2,b)*VxNx(0,a) - - r23*NxFi(2,a)*VxNx(0,b)))*mu*Jac; - Ku = w*af*(T1 + T2 + T3 + Tv + BtDB); + Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(6,a,b)); lKd(6,a,b) = lKd(6,a,b) + Ku; T2 = af*Jac*tauC*rho*NxFi(2,a)*NxFi(0,b); - Tv = af*mu*Jac*(NxFi(0,a)*NxFi(2,b) - - r23*NxFi(2,a)*NxFi(0,b)); + Tv = af*Kvis_v(6,a,b); lK(8,a,b) = lK(8,a,b) + w*(T2 + Tv) + afm*Ku; // dM_3/dV_2 + af/am *dM_3/dU_2 @@ -1607,17 +1524,12 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T1 = Jac*rho*vd(2)*Nw(a)*NxFi(1,b); T2 = -tauC*Jac*NxFi(2,a)*VxNx(1,b); T3 = Jac*rCl*(NxFi(2,a)*NxFi(1,b) - NxFi(1,a)*NxFi(2,b)); - - Tv = (2.0*(DdNx(2,a)*NxFi(1,b) - DdNx(2,b)*NxFi(1,a)) - - (NxNx*VxFi(2,1) + NxFi(2,b)*VxNx(1,a) - - r23*NxFi(2,a)*VxNx(1,b)))*mu*Jac; - Ku = w*af*(T1 + T2 + T3 + Tv + BtDB); + Ku = w*af*(T1 + T2 + T3 + BtDB + Kvis_u(7,a,b)); lKd(7,a,b) = lKd(7,a,b) + Ku; T2 = af*Jac*tauC*rho*NxFi(2,a)*NxFi(1,b); - Tv = af*mu*Jac*(NxFi(1,a)*NxFi(2,b) - - r23*NxFi(2,a)*NxFi(1,b)); + Tv = af*Kvis_v(7,a,b); lK(9,a,b) = lK(9,a,b) + w*(T2 + Tv) + afm*Ku; @@ -1629,17 +1541,13 @@ void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in T1 = Jac*rho*vd(2)*Nw(a)*NxFi(2,b); T2 = -tauC*Jac*NxFi(2,a)*VxNx(2,b); - - Tv = (2.0*(DdNx(2,a)*NxFi(2,b) - DdNx(2,b)*NxFi(2,a)) - - (NxNx*VxFi(2,2) + NxFi(2,b)*VxNx(2,a) - - r23*NxFi(2,a)*VxNx(2,b)))*mu*Jac; - Ku = w*af*(T1 + T2 + Tv + BtDB + NxSNx); + Ku = w*af*(T1 + T2 + BtDB + NxSNx + Kvis_u(8,a,b)); lKd(8,a,b) = lKd(8,a,b) + Ku; T1 = am*Jac*rho*Nw(a)*Nw(b); T2 = T1 + af*Jac*tauC*rho*NxFi(2,a)*NxFi(2,b); - Tv = af*mu*Jac*(r13*NxFi(2,a)*NxFi(2,b) + NxNx); + Tv = af*Kvis_v(8,a,b); lK(10,a,b) = lK(10,a,b) + w*(T2 + Tv) + afm*Ku; } diff --git a/Code/Source/svFSI/ustruct.h b/Code/Source/svFSI/ustruct.h index 972e88a0..71698323 100644 --- a/Code/Source/svFSI/ustruct.h +++ b/Code/Source/svFSI/ustruct.h @@ -51,7 +51,7 @@ int get_col_ptr(ComMod& com_mod, const int rowN, const int colN); void ustruct_2d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, const double w, const double Je, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& al, const Array& yl, - const Array& dl, const Array& bfl, const Array& Kxi, Array& lR, Array3& lK, + const Array& dl, const Array& bfl, Array& lR, Array3& lK, Array3& lKd); void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, @@ -63,7 +63,7 @@ void ustruct_2d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const in void ustruct_3d_c(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, const double w, const double Je, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& al, const Array& yl, - const Array& dl, const Array& bfl, const Array& Kxi, Array& lR, Array3& lK, + const Array& dl, const Array& bfl, Array& lR, Array3& lK, Array3& lKd); void ustruct_3d_m(ComMod& com_mod, CepMod& cep_mod, const bool vmsFlag, const int eNoNw, const int eNoNq, diff --git a/tests/cases/stokes/manufactured_solution/P1P1/N008/svFSIplus.xml b/tests/cases/stokes/manufactured_solution/P1P1/N008/svFSIplus.xml index 9f7c24e5..4f0e2837 100644 --- a/tests/cases/stokes/manufactured_solution/P1P1/N008/svFSIplus.xml +++ b/tests/cases/stokes/manufactured_solution/P1P1/N008/svFSIplus.xml @@ -47,7 +47,7 @@ 1e-12 false - + "Constant" > 1.0 diff --git a/tests/cases/stokes/manufactured_solution/P2P1/N016/svFSIplus.xml b/tests/cases/stokes/manufactured_solution/P2P1/N016/svFSIplus.xml index 54288f9b..bbc80b4b 100644 --- a/tests/cases/stokes/manufactured_solution/P2P1/N016/svFSIplus.xml +++ b/tests/cases/stokes/manufactured_solution/P2P1/N016/svFSIplus.xml @@ -47,7 +47,7 @@ 1e-9 true - + "Constant" > 1.0 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/README.md b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/README.md new file mode 100644 index 00000000..89e11caa --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/README.md @@ -0,0 +1,32 @@ +This test case simulates pulling and releasing a slab of material +described by the Guccione material model and a Newtonian solid viscosity model. +In this viscosity model, the viscous deviatoric Cauchy stress is identical to that +for a Newtonian fluid: + +$$ +\sigma^{dev}_{vis} = 2 \mu \mathbf{d}^{dev} +$$ +where +$$ +\mathbf{d}^{dev} = \frac{1}{2} (\nabla_x \mathbf{v} + (\nabla_x \mathbf{v})^T) - \frac{1}{3} (\nabla_x \cdot \mathbf{v}) \mathbf{I} +$$ + +The viscous part of the 2nd Piola-Kirchhoff stress is then given by a pull-back operation +$$ +\mathbf{S}_{vis} = 2 \mu J \mathbf{F}^{-1} \mathbf{d}^{dev} \mathbf{F}^{-T} +$$ + +The load profile is a 0.25s ramp to the max value, then held for 0.25s, then a 0.25s ramp down to 0, then held for 0.25s. +The load is applied on the Z1 face in the z-direction. + +![Load Profile](load.png) + +The load data is defined in `load.dat`, which can be generated with +`load.py` + +The resulting deformation is shown in the video below: +![Deformation](movie.gif) + +The z-displacement of the Z1 face is plotted below: +![Z1 Displacement](z1_face_displacement.png) + diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/generate_load.py b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/generate_load.py new file mode 100644 index 00000000..88e2ac57 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/generate_load.py @@ -0,0 +1,34 @@ +import numpy as np +import os + +# Go to directory of this script +os.chdir(os.path.dirname(os.path.realpath(__file__))) + +# Number of timesteps and number of Fourier modes +n_timesteps = 1001 +n_modes = 256 + +# Generate time values from 0 to 1 +time = np.linspace(0, 1, n_timesteps) + +# Generate 0.25s ramp to max load, then hold for 0.25s, then 0.25s ramp down to 0, then hold for 0.25s +load = np.zeros(n_timesteps) +max_load= -1e2 +load[time < 0.25] = max_load * time[time < 0.25] / 0.25 +load[(time >= 0.25) & (time < 0.5)] = max_load +load[(time >= 0.5) & (time < 0.75)] = max_load * (1 - (time[(time >= 0.5) & (time < 0.75)] - 0.5) / 0.25) +load[time >= 0.75] = 0 + +# Write the time and load values to a text file +with open("load.dat", "w") as file: + file.write(f"{n_timesteps} {n_modes}\n") + for t, s in zip(time, load): + file.write(f"{t:.3f} {s:.3f}\n") + +# Plot the load values +import matplotlib.pyplot as plt +plt.plot(time, load) +plt.xlabel("Time (s)") +plt.ylabel("load (dynes/cm^2)") +plt.title("Active load") +plt.savefig("load.png") \ No newline at end of file diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/load.dat b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/load.dat new file mode 100644 index 00000000..0e638709 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/load.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1baa9ed6cf748d42b962453326733347bb3f82c1bd9c49a653568896e15e1e6c +size 13723 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/load.png b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/load.png new file mode 100644 index 00000000..5b70832a --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/load.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c58557b9830a4b793a2bb89cf130bdc4540104d8190db0c64dc59e0da1cf8d9b +size 22901 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-complete.mesh.vtu b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-complete.mesh.vtu new file mode 100644 index 00000000..f09242eb --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d74418f571657908b5c4e73f7cc5ce040fd300ac444b1c893ab23001bb6e09ae +size 7499 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/X0.vtp b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/X0.vtp new file mode 100644 index 00000000..904401ad --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/X0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3f71480fecee6f177719ab68b14bb2b61bbbce8fc745c87576fd51f3fb0b3ef8 +size 3655 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/X1.vtp b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/X1.vtp new file mode 100644 index 00000000..20b2d0e5 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/X1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bd46bcf90705c0e062aa658846aacb9091334a56fffac81715631628cd1f097f +size 3669 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Y0.vtp b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Y0.vtp new file mode 100644 index 00000000..19c7f008 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Y0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1ef44a9a441563692f58baccf576b3c10e6df4277046d1054db5ab821d498bf0 +size 4795 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Y1.vtp b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Y1.vtp new file mode 100644 index 00000000..7ac40cb0 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Y1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4d2ad3c14f29583d9d2015b0c0ba681b66c0d2380e4b8851ca30a2cb81773183 +size 4629 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Z0.vtp b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Z0.vtp new file mode 100644 index 00000000..49112c41 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Z0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6ee78716770c3ed486e5e81d32f90226e215bc095b5cda4f5c708218b84ee407 +size 3188 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Z1.vtp b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Z1.vtp new file mode 100644 index 00000000..22122b6f --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Z1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c2bdc7813853e641c5a71c20f1cea215c9c0215523f51272c5650077614dcb29 +size 3245 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/movie.avi b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/movie.avi new file mode 100644 index 00000000..cc571072 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/movie.avi @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bb7fff34151f35d1e773ffa24943fc829e66b759ea05c0e7e06862fb36c297b0 +size 6689058 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/movie.gif b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/movie.gif new file mode 100644 index 00000000..bcb06b5a --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/movie.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8248b3481ce5c4614ab8c915c370b2d0f479b11a4ea4d6964c8ed09f7e726935 +size 915192 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/plot_z1_face_displacement.py b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/plot_z1_face_displacement.py new file mode 100644 index 00000000..609cb008 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/plot_z1_face_displacement.py @@ -0,0 +1,145 @@ +''' +Plots the displacement of the z1 face of the mesh over time +''' + +import matplotlib.pyplot as plt +import numpy as np +import pyvista as pv +import os +import glob +import xml.etree.ElementTree as ET + +def get_timestep(result_file): + ''' + Extracts the timestep of a results_###.vtu file + ''' + # If file path is provided, get the file name + file_name = os.path.basename(result_file) + + # Get the ###.vtu part + s = file_name.split('_')[-1] + + # Get the ### part + s = s.split('.')[0] + + # Return as integer + return int(s) + +def get_start_end_step(results_folder): + """ + Automatically determine the start timestep, end timestep, and step size based on all + svFSI results file in results_folder. + + Args: + results_folder: A string of absolute file path of folder containing results of + svFSIplus simulation. This usually ends with 16-procs/ or other + number. + + Returns: + (start_timestep, end_timestep, step): A tuple of 3 integers, giving the first + timestep of results to process, the last timestep, and the step size. + + """ + + # Get list of all .vtu files in results_folder sorted by time step + list_of_files = sorted( filter( os.path.isfile, + glob.glob(os.path.join(results_folder, '*.vtu')) ), key = get_timestep) + + # Get start time from the first result file (list_of_files[0]) + start_file_name = os.path.basename(list_of_files[0]) + start_timestep = int("".join([i for i in start_file_name if i.isdigit()])) + print('Start timestep:', start_timestep) + + # Get end time from the last result file (list_of_files[-1]) + end_file_name = os.path.basename(list_of_files[-1]) + end_timestep = int("".join([i for i in end_file_name if i.isdigit()])) + print('End timestep:', end_timestep) + + # Get step size by looking at second time step + start_plus_one_file_name = os.path.basename(list_of_files[1]) + start_time_plus_one = int("".join([i for i in start_plus_one_file_name if i.isdigit()])) + step = start_time_plus_one - start_timestep + print('Step:', step) + + return (start_timestep, end_timestep, step) + +def get_timestep_size(svfsiplus_xml_file): + ''' + Reads the timestep size from an svFSIplus XML file. + + Args: + svfsiplus_xml_file: The file path to the svFSIplus XML input file + + Returns: + timestep_size: The size of the timestep in the input file + ''' + # Parse the XML file + tree = ET.parse(svfsiplus_xml_file) + root = tree.getroot() + + # Find the element containing the Time_step_size + timestep_size = None + for elem in root.iter(): + if elem.tag == 'Time_step_size': + timestep_size = float(elem.text) + break + + return timestep_size + + +# Move to this script's directory +os.chdir(os.path.dirname(os.path.realpath(__file__))) + +# Load the z1 face +mesh = pv.read("mesh/mesh-surfaces/Z1.vtp") + +# Get the start timestep, end timestep, and step size +start_timestep, end_timestep, step = get_start_end_step("4-procs/") + +# Get the timestep size +timestep_size = get_timestep_size("svFSIPlus.xml") + +# Get the time values +time = np.arange(start_timestep, end_timestep + 1, step) * timestep_size + +# Loop over the results files +displacements = [] +for n in range(start_timestep, end_timestep + 1, step): + # Load the results file + result_file = f"4-procs/result_{n:03d}.vtu" + result = pv.read(result_file) + + # Sample displacement at the z1 face + displacement = mesh.sample(result).point_data['Displacement'][:, 2].mean() + + # Append to the list + displacements.append(displacement) + +# Save the displacements to a text file +with open("z1_face_displacement.dat", "w") as file: + file.write(f"{len(time)}\n") + for t, s in zip(time, displacements): + file.write(f"{t:.3f} {s:.3f}\n") + + +# Plot the displacement +plt.plot(time, displacements, label="With viscosity") +plt.xlabel("Time (s)") +plt.ylabel("z Displacement (cm)") +plt.title("Z1 face displacement") + +# Plot the displacement with no viscosity +if os.path.exists("z1_face_displacement_no_viscosity.dat"): + with open("z1_face_displacement_no_viscosity.dat", "r") as file: + lines = file.readlines() + n = int(lines[0]) + time_no_viscosity = np.zeros(n) + displacements_no_viscosity = np.zeros(n) + for i, line in enumerate(lines[1:]): + time_no_viscosity[i], displacements_no_viscosity[i] = map(float, line.split()) + + plt.plot(time_no_viscosity, displacements_no_viscosity, label="Without viscosity") + + +plt.legend() +plt.savefig("z1_face_displacement.png") \ No newline at end of file diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/result_001.vtu b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/result_001.vtu new file mode 100644 index 00000000..05be8a51 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/result_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:88da8f9d0b6275eef0731efd4774edb3dd32acc2e8f3cde954c2c7b4f87e4c90 +size 76626 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/svFSIplus.xml b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/svFSIplus.xml new file mode 100644 index 00000000..8784bd6d --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/svFSIplus.xml @@ -0,0 +1,125 @@ + + + + + 0 + 3 + 1 + 0.01 + 0.50 + STOP_SIM + + 1 + result + 1 + 1 + + 100 + 0 + + 1 + 0 + 0 + + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/X0.vtp + + + + mesh/mesh-surfaces/X1.vtp + + + + mesh/mesh-surfaces/Y0.vtp + + + + mesh/mesh-surfaces/Y1.vtp + + + + mesh/mesh-surfaces/Z0.vtp + + + + mesh/mesh-surfaces/Z1.vtp + + + (0.0, 0.0, 1.0) + (1.0, 0.0, 0.0) + + + + + + + + + + true + 1 + 12 + 1e-9 + + 1.0e-2 + 1.0e5 + 0.48333 + + + 50.0 + + + + 100.0 + 8.0 + 6.0 + 12.0 + + + ST91 + + + true + true + true + true + true + true + true + true + + + + + fsils + + 1e-10 + 400 + + + + Dir + 0.0 + true + false + + + + Neu + Unsteady + load.dat + false + + + + + + + diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement.dat b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement.dat new file mode 100644 index 00000000..8536ab01 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8b9be491d17fad5eac0f4f834a8e9f26be1d13f70937fd6b83ae34273ea4260a +size 1204 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement.png b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement.png new file mode 100644 index 00000000..1d22a601 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:23fa1813c8c831ad04236a3ccaedce2119e510a3b99dd67e06b13a5f1dbc18c5 +size 45469 diff --git a/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement_no_viscosity.dat b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement_no_viscosity.dat new file mode 100644 index 00000000..3ef76690 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement_no_viscosity.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2c29c7f8626c2e0132ded468a36acc9c19f1a3142c49aac9168795d3398739ac +size 1219 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/README.md b/tests/cases/struct/tensile_adventitia_Potential_viscosity/README.md new file mode 100644 index 00000000..b7701399 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/README.md @@ -0,0 +1,28 @@ +This test case simulates pulling and releasing a slab of material +described by the Guccione material model and a pseudo-potential solid viscosity model. +In this viscosity model, the so-called viscous pseudo-potential $\Psi_{vis}$ is +given by +$$ +\Psi_{vis} = \frac{\mu}{2} \text{tr}(\dot{\mathbf{E}}) +$$ +The viscous part of the 2nd Piola-Kirchhoff stress is then given by +$$ +S_{vis} = \frac{\partial \Psi_{vis}}{\partial \dot{\mathbf{E}}} = \mu \dot{\mathbf{E}} +$$ + + + +The load profile is a 0.25s ramp to the max value, then held for 0.25s, then a 0.25s ramp down to 0, then held for 0.25s. +The load is applied on the Z1 face in the z-direction. + +![Load Profile](load.png) + +The load data is defined in `load.dat`, which can be generated with +`load.py` + +The resulting deformation is shown in the video below: +![Deformation](movie.gif) + +The z-displacement of the Z1 face is plotted below: +![Z1 Displacement](z1_face_displacement.png) + diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/generate_load.py b/tests/cases/struct/tensile_adventitia_Potential_viscosity/generate_load.py new file mode 100644 index 00000000..88e2ac57 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/generate_load.py @@ -0,0 +1,34 @@ +import numpy as np +import os + +# Go to directory of this script +os.chdir(os.path.dirname(os.path.realpath(__file__))) + +# Number of timesteps and number of Fourier modes +n_timesteps = 1001 +n_modes = 256 + +# Generate time values from 0 to 1 +time = np.linspace(0, 1, n_timesteps) + +# Generate 0.25s ramp to max load, then hold for 0.25s, then 0.25s ramp down to 0, then hold for 0.25s +load = np.zeros(n_timesteps) +max_load= -1e2 +load[time < 0.25] = max_load * time[time < 0.25] / 0.25 +load[(time >= 0.25) & (time < 0.5)] = max_load +load[(time >= 0.5) & (time < 0.75)] = max_load * (1 - (time[(time >= 0.5) & (time < 0.75)] - 0.5) / 0.25) +load[time >= 0.75] = 0 + +# Write the time and load values to a text file +with open("load.dat", "w") as file: + file.write(f"{n_timesteps} {n_modes}\n") + for t, s in zip(time, load): + file.write(f"{t:.3f} {s:.3f}\n") + +# Plot the load values +import matplotlib.pyplot as plt +plt.plot(time, load) +plt.xlabel("Time (s)") +plt.ylabel("load (dynes/cm^2)") +plt.title("Active load") +plt.savefig("load.png") \ No newline at end of file diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/load.dat b/tests/cases/struct/tensile_adventitia_Potential_viscosity/load.dat new file mode 100644 index 00000000..0e638709 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/load.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1baa9ed6cf748d42b962453326733347bb3f82c1bd9c49a653568896e15e1e6c +size 13723 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/load.png b/tests/cases/struct/tensile_adventitia_Potential_viscosity/load.png new file mode 100644 index 00000000..5b70832a --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/load.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c58557b9830a4b793a2bb89cf130bdc4540104d8190db0c64dc59e0da1cf8d9b +size 22901 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-complete.mesh.vtu b/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-complete.mesh.vtu new file mode 100644 index 00000000..f09242eb --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d74418f571657908b5c4e73f7cc5ce040fd300ac444b1c893ab23001bb6e09ae +size 7499 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/X0.vtp b/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/X0.vtp new file mode 100644 index 00000000..904401ad --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/X0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3f71480fecee6f177719ab68b14bb2b61bbbce8fc745c87576fd51f3fb0b3ef8 +size 3655 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/X1.vtp b/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/X1.vtp new file mode 100644 index 00000000..20b2d0e5 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/X1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bd46bcf90705c0e062aa658846aacb9091334a56fffac81715631628cd1f097f +size 3669 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Y0.vtp b/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Y0.vtp new file mode 100644 index 00000000..19c7f008 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Y0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1ef44a9a441563692f58baccf576b3c10e6df4277046d1054db5ab821d498bf0 +size 4795 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Y1.vtp b/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Y1.vtp new file mode 100644 index 00000000..7ac40cb0 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Y1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4d2ad3c14f29583d9d2015b0c0ba681b66c0d2380e4b8851ca30a2cb81773183 +size 4629 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Z0.vtp b/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Z0.vtp new file mode 100644 index 00000000..49112c41 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Z0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6ee78716770c3ed486e5e81d32f90226e215bc095b5cda4f5c708218b84ee407 +size 3188 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Z1.vtp b/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Z1.vtp new file mode 100644 index 00000000..22122b6f --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Z1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c2bdc7813853e641c5a71c20f1cea215c9c0215523f51272c5650077614dcb29 +size 3245 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/movie.avi b/tests/cases/struct/tensile_adventitia_Potential_viscosity/movie.avi new file mode 100644 index 00000000..31bbc25c --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/movie.avi @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ba4cd361b458769225571014243263bcfc0da4cd823f33f127504b6be3277c76 +size 6617286 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/movie.gif b/tests/cases/struct/tensile_adventitia_Potential_viscosity/movie.gif new file mode 100644 index 00000000..f5cfd62a --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/movie.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:cd29195603702282ef78ddb057f52f9deebebe68d2e7333e218d8744aaac5236 +size 4112863 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/plot_z1_face_displacement.py b/tests/cases/struct/tensile_adventitia_Potential_viscosity/plot_z1_face_displacement.py new file mode 100644 index 00000000..609cb008 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/plot_z1_face_displacement.py @@ -0,0 +1,145 @@ +''' +Plots the displacement of the z1 face of the mesh over time +''' + +import matplotlib.pyplot as plt +import numpy as np +import pyvista as pv +import os +import glob +import xml.etree.ElementTree as ET + +def get_timestep(result_file): + ''' + Extracts the timestep of a results_###.vtu file + ''' + # If file path is provided, get the file name + file_name = os.path.basename(result_file) + + # Get the ###.vtu part + s = file_name.split('_')[-1] + + # Get the ### part + s = s.split('.')[0] + + # Return as integer + return int(s) + +def get_start_end_step(results_folder): + """ + Automatically determine the start timestep, end timestep, and step size based on all + svFSI results file in results_folder. + + Args: + results_folder: A string of absolute file path of folder containing results of + svFSIplus simulation. This usually ends with 16-procs/ or other + number. + + Returns: + (start_timestep, end_timestep, step): A tuple of 3 integers, giving the first + timestep of results to process, the last timestep, and the step size. + + """ + + # Get list of all .vtu files in results_folder sorted by time step + list_of_files = sorted( filter( os.path.isfile, + glob.glob(os.path.join(results_folder, '*.vtu')) ), key = get_timestep) + + # Get start time from the first result file (list_of_files[0]) + start_file_name = os.path.basename(list_of_files[0]) + start_timestep = int("".join([i for i in start_file_name if i.isdigit()])) + print('Start timestep:', start_timestep) + + # Get end time from the last result file (list_of_files[-1]) + end_file_name = os.path.basename(list_of_files[-1]) + end_timestep = int("".join([i for i in end_file_name if i.isdigit()])) + print('End timestep:', end_timestep) + + # Get step size by looking at second time step + start_plus_one_file_name = os.path.basename(list_of_files[1]) + start_time_plus_one = int("".join([i for i in start_plus_one_file_name if i.isdigit()])) + step = start_time_plus_one - start_timestep + print('Step:', step) + + return (start_timestep, end_timestep, step) + +def get_timestep_size(svfsiplus_xml_file): + ''' + Reads the timestep size from an svFSIplus XML file. + + Args: + svfsiplus_xml_file: The file path to the svFSIplus XML input file + + Returns: + timestep_size: The size of the timestep in the input file + ''' + # Parse the XML file + tree = ET.parse(svfsiplus_xml_file) + root = tree.getroot() + + # Find the element containing the Time_step_size + timestep_size = None + for elem in root.iter(): + if elem.tag == 'Time_step_size': + timestep_size = float(elem.text) + break + + return timestep_size + + +# Move to this script's directory +os.chdir(os.path.dirname(os.path.realpath(__file__))) + +# Load the z1 face +mesh = pv.read("mesh/mesh-surfaces/Z1.vtp") + +# Get the start timestep, end timestep, and step size +start_timestep, end_timestep, step = get_start_end_step("4-procs/") + +# Get the timestep size +timestep_size = get_timestep_size("svFSIPlus.xml") + +# Get the time values +time = np.arange(start_timestep, end_timestep + 1, step) * timestep_size + +# Loop over the results files +displacements = [] +for n in range(start_timestep, end_timestep + 1, step): + # Load the results file + result_file = f"4-procs/result_{n:03d}.vtu" + result = pv.read(result_file) + + # Sample displacement at the z1 face + displacement = mesh.sample(result).point_data['Displacement'][:, 2].mean() + + # Append to the list + displacements.append(displacement) + +# Save the displacements to a text file +with open("z1_face_displacement.dat", "w") as file: + file.write(f"{len(time)}\n") + for t, s in zip(time, displacements): + file.write(f"{t:.3f} {s:.3f}\n") + + +# Plot the displacement +plt.plot(time, displacements, label="With viscosity") +plt.xlabel("Time (s)") +plt.ylabel("z Displacement (cm)") +plt.title("Z1 face displacement") + +# Plot the displacement with no viscosity +if os.path.exists("z1_face_displacement_no_viscosity.dat"): + with open("z1_face_displacement_no_viscosity.dat", "r") as file: + lines = file.readlines() + n = int(lines[0]) + time_no_viscosity = np.zeros(n) + displacements_no_viscosity = np.zeros(n) + for i, line in enumerate(lines[1:]): + time_no_viscosity[i], displacements_no_viscosity[i] = map(float, line.split()) + + plt.plot(time_no_viscosity, displacements_no_viscosity, label="Without viscosity") + + +plt.legend() +plt.savefig("z1_face_displacement.png") \ No newline at end of file diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/result_001.vtu b/tests/cases/struct/tensile_adventitia_Potential_viscosity/result_001.vtu new file mode 100644 index 00000000..d18a1ea6 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/result_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:238e6783fbe88d0453f5ea64cd615c7eff9225a63a545d468844f339895c31f2 +size 77010 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/svFSIplus.xml b/tests/cases/struct/tensile_adventitia_Potential_viscosity/svFSIplus.xml new file mode 100644 index 00000000..9b1eeeaa --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/svFSIplus.xml @@ -0,0 +1,125 @@ + + + + + 0 + 3 + 1 + 0.01 + 0.50 + STOP_SIM + + 1 + result + 1 + 1 + + 100 + 0 + + 1 + 0 + 0 + + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/X0.vtp + + + + mesh/mesh-surfaces/X1.vtp + + + + mesh/mesh-surfaces/Y0.vtp + + + + mesh/mesh-surfaces/Y1.vtp + + + + mesh/mesh-surfaces/Z0.vtp + + + + mesh/mesh-surfaces/Z1.vtp + + + (0.0, 0.0, 1.0) + (1.0, 0.0, 0.0) + + + + + + + + + + true + 1 + 12 + 1e-9 + + 1.0e-2 + 1.0e5 + 0.48333 + + + 50.0 + + + + 100.0 + 8.0 + 6.0 + 12.0 + + + ST91 + + + true + true + true + true + true + true + true + true + + + + + fsils + + 1e-10 + 400 + + + + Dir + 0.0 + true + false + + + + Neu + Unsteady + load.dat + false + + + + + + + diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/z1_face_displacement.dat b/tests/cases/struct/tensile_adventitia_Potential_viscosity/z1_face_displacement.dat new file mode 100644 index 00000000..6bfbd991 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/z1_face_displacement.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:edd0644b1495f7143f126d05761955e683c52b3629c7ecdbfe145536b2cdd496 +size 1204 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/z1_face_displacement.png b/tests/cases/struct/tensile_adventitia_Potential_viscosity/z1_face_displacement.png new file mode 100644 index 00000000..fe2f4ae2 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/z1_face_displacement.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6619455ecc4d4ea2dbcf3059ed3d547e13a41520d89ecd2a1dc7ed194ba2adc2 +size 44793 diff --git a/tests/cases/struct/tensile_adventitia_Potential_viscosity/z1_face_displacement_no_viscosity.dat b/tests/cases/struct/tensile_adventitia_Potential_viscosity/z1_face_displacement_no_viscosity.dat new file mode 100644 index 00000000..3ef76690 --- /dev/null +++ b/tests/cases/struct/tensile_adventitia_Potential_viscosity/z1_face_displacement_no_viscosity.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2c29c7f8626c2e0132ded468a36acc9c19f1a3142c49aac9168795d3398739ac +size 1219 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/README.md b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/README.md new file mode 100644 index 00000000..89e11caa --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/README.md @@ -0,0 +1,32 @@ +This test case simulates pulling and releasing a slab of material +described by the Guccione material model and a Newtonian solid viscosity model. +In this viscosity model, the viscous deviatoric Cauchy stress is identical to that +for a Newtonian fluid: + +$$ +\sigma^{dev}_{vis} = 2 \mu \mathbf{d}^{dev} +$$ +where +$$ +\mathbf{d}^{dev} = \frac{1}{2} (\nabla_x \mathbf{v} + (\nabla_x \mathbf{v})^T) - \frac{1}{3} (\nabla_x \cdot \mathbf{v}) \mathbf{I} +$$ + +The viscous part of the 2nd Piola-Kirchhoff stress is then given by a pull-back operation +$$ +\mathbf{S}_{vis} = 2 \mu J \mathbf{F}^{-1} \mathbf{d}^{dev} \mathbf{F}^{-T} +$$ + +The load profile is a 0.25s ramp to the max value, then held for 0.25s, then a 0.25s ramp down to 0, then held for 0.25s. +The load is applied on the Z1 face in the z-direction. + +![Load Profile](load.png) + +The load data is defined in `load.dat`, which can be generated with +`load.py` + +The resulting deformation is shown in the video below: +![Deformation](movie.gif) + +The z-displacement of the Z1 face is plotted below: +![Z1 Displacement](z1_face_displacement.png) + diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/generate_load.py b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/generate_load.py new file mode 100644 index 00000000..88e2ac57 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/generate_load.py @@ -0,0 +1,34 @@ +import numpy as np +import os + +# Go to directory of this script +os.chdir(os.path.dirname(os.path.realpath(__file__))) + +# Number of timesteps and number of Fourier modes +n_timesteps = 1001 +n_modes = 256 + +# Generate time values from 0 to 1 +time = np.linspace(0, 1, n_timesteps) + +# Generate 0.25s ramp to max load, then hold for 0.25s, then 0.25s ramp down to 0, then hold for 0.25s +load = np.zeros(n_timesteps) +max_load= -1e2 +load[time < 0.25] = max_load * time[time < 0.25] / 0.25 +load[(time >= 0.25) & (time < 0.5)] = max_load +load[(time >= 0.5) & (time < 0.75)] = max_load * (1 - (time[(time >= 0.5) & (time < 0.75)] - 0.5) / 0.25) +load[time >= 0.75] = 0 + +# Write the time and load values to a text file +with open("load.dat", "w") as file: + file.write(f"{n_timesteps} {n_modes}\n") + for t, s in zip(time, load): + file.write(f"{t:.3f} {s:.3f}\n") + +# Plot the load values +import matplotlib.pyplot as plt +plt.plot(time, load) +plt.xlabel("Time (s)") +plt.ylabel("load (dynes/cm^2)") +plt.title("Active load") +plt.savefig("load.png") \ No newline at end of file diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/load.dat b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/load.dat new file mode 100644 index 00000000..0e638709 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/load.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1baa9ed6cf748d42b962453326733347bb3f82c1bd9c49a653568896e15e1e6c +size 13723 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/load.png b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/load.png new file mode 100644 index 00000000..5b70832a --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/load.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c58557b9830a4b793a2bb89cf130bdc4540104d8190db0c64dc59e0da1cf8d9b +size 22901 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-complete.mesh.vtu b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-complete.mesh.vtu new file mode 100644 index 00000000..f09242eb --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d74418f571657908b5c4e73f7cc5ce040fd300ac444b1c893ab23001bb6e09ae +size 7499 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/X0.vtp b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/X0.vtp new file mode 100644 index 00000000..904401ad --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/X0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3f71480fecee6f177719ab68b14bb2b61bbbce8fc745c87576fd51f3fb0b3ef8 +size 3655 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/X1.vtp b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/X1.vtp new file mode 100644 index 00000000..20b2d0e5 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/X1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bd46bcf90705c0e062aa658846aacb9091334a56fffac81715631628cd1f097f +size 3669 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Y0.vtp b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Y0.vtp new file mode 100644 index 00000000..19c7f008 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Y0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1ef44a9a441563692f58baccf576b3c10e6df4277046d1054db5ab821d498bf0 +size 4795 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Y1.vtp b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Y1.vtp new file mode 100644 index 00000000..7ac40cb0 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Y1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4d2ad3c14f29583d9d2015b0c0ba681b66c0d2380e4b8851ca30a2cb81773183 +size 4629 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Z0.vtp b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Z0.vtp new file mode 100644 index 00000000..49112c41 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Z0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6ee78716770c3ed486e5e81d32f90226e215bc095b5cda4f5c708218b84ee407 +size 3188 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Z1.vtp b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Z1.vtp new file mode 100644 index 00000000..22122b6f --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/mesh/mesh-surfaces/Z1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c2bdc7813853e641c5a71c20f1cea215c9c0215523f51272c5650077614dcb29 +size 3245 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/movie.avi b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/movie.avi new file mode 100644 index 00000000..7aa01e56 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/movie.avi @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6d3b3aeaa3f7ed69a0380f2ef521b469845948d0c73a40ddbfe2c525cc7d6796 +size 6696188 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/movie.gif b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/movie.gif new file mode 100644 index 00000000..abf13afb --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/movie.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:70ec054e059b6f00ab53239a4230f9c77934aa25583144aeb017d9644993265e +size 938216 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/plot_z1_face_displacement.py b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/plot_z1_face_displacement.py new file mode 100644 index 00000000..609cb008 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/plot_z1_face_displacement.py @@ -0,0 +1,145 @@ +''' +Plots the displacement of the z1 face of the mesh over time +''' + +import matplotlib.pyplot as plt +import numpy as np +import pyvista as pv +import os +import glob +import xml.etree.ElementTree as ET + +def get_timestep(result_file): + ''' + Extracts the timestep of a results_###.vtu file + ''' + # If file path is provided, get the file name + file_name = os.path.basename(result_file) + + # Get the ###.vtu part + s = file_name.split('_')[-1] + + # Get the ### part + s = s.split('.')[0] + + # Return as integer + return int(s) + +def get_start_end_step(results_folder): + """ + Automatically determine the start timestep, end timestep, and step size based on all + svFSI results file in results_folder. + + Args: + results_folder: A string of absolute file path of folder containing results of + svFSIplus simulation. This usually ends with 16-procs/ or other + number. + + Returns: + (start_timestep, end_timestep, step): A tuple of 3 integers, giving the first + timestep of results to process, the last timestep, and the step size. + + """ + + # Get list of all .vtu files in results_folder sorted by time step + list_of_files = sorted( filter( os.path.isfile, + glob.glob(os.path.join(results_folder, '*.vtu')) ), key = get_timestep) + + # Get start time from the first result file (list_of_files[0]) + start_file_name = os.path.basename(list_of_files[0]) + start_timestep = int("".join([i for i in start_file_name if i.isdigit()])) + print('Start timestep:', start_timestep) + + # Get end time from the last result file (list_of_files[-1]) + end_file_name = os.path.basename(list_of_files[-1]) + end_timestep = int("".join([i for i in end_file_name if i.isdigit()])) + print('End timestep:', end_timestep) + + # Get step size by looking at second time step + start_plus_one_file_name = os.path.basename(list_of_files[1]) + start_time_plus_one = int("".join([i for i in start_plus_one_file_name if i.isdigit()])) + step = start_time_plus_one - start_timestep + print('Step:', step) + + return (start_timestep, end_timestep, step) + +def get_timestep_size(svfsiplus_xml_file): + ''' + Reads the timestep size from an svFSIplus XML file. + + Args: + svfsiplus_xml_file: The file path to the svFSIplus XML input file + + Returns: + timestep_size: The size of the timestep in the input file + ''' + # Parse the XML file + tree = ET.parse(svfsiplus_xml_file) + root = tree.getroot() + + # Find the element containing the Time_step_size + timestep_size = None + for elem in root.iter(): + if elem.tag == 'Time_step_size': + timestep_size = float(elem.text) + break + + return timestep_size + + +# Move to this script's directory +os.chdir(os.path.dirname(os.path.realpath(__file__))) + +# Load the z1 face +mesh = pv.read("mesh/mesh-surfaces/Z1.vtp") + +# Get the start timestep, end timestep, and step size +start_timestep, end_timestep, step = get_start_end_step("4-procs/") + +# Get the timestep size +timestep_size = get_timestep_size("svFSIPlus.xml") + +# Get the time values +time = np.arange(start_timestep, end_timestep + 1, step) * timestep_size + +# Loop over the results files +displacements = [] +for n in range(start_timestep, end_timestep + 1, step): + # Load the results file + result_file = f"4-procs/result_{n:03d}.vtu" + result = pv.read(result_file) + + # Sample displacement at the z1 face + displacement = mesh.sample(result).point_data['Displacement'][:, 2].mean() + + # Append to the list + displacements.append(displacement) + +# Save the displacements to a text file +with open("z1_face_displacement.dat", "w") as file: + file.write(f"{len(time)}\n") + for t, s in zip(time, displacements): + file.write(f"{t:.3f} {s:.3f}\n") + + +# Plot the displacement +plt.plot(time, displacements, label="With viscosity") +plt.xlabel("Time (s)") +plt.ylabel("z Displacement (cm)") +plt.title("Z1 face displacement") + +# Plot the displacement with no viscosity +if os.path.exists("z1_face_displacement_no_viscosity.dat"): + with open("z1_face_displacement_no_viscosity.dat", "r") as file: + lines = file.readlines() + n = int(lines[0]) + time_no_viscosity = np.zeros(n) + displacements_no_viscosity = np.zeros(n) + for i, line in enumerate(lines[1:]): + time_no_viscosity[i], displacements_no_viscosity[i] = map(float, line.split()) + + plt.plot(time_no_viscosity, displacements_no_viscosity, label="Without viscosity") + + +plt.legend() +plt.savefig("z1_face_displacement.png") \ No newline at end of file diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/result_001.vtu b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/result_001.vtu new file mode 100644 index 00000000..13b5d36f --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/result_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:56a59efe0cb8214f54659f34cac97ef0c5028652abea1f434d23037768b9a26f +size 81000 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/svFSIplus.xml b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/svFSIplus.xml new file mode 100644 index 00000000..33b8c509 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/svFSIplus.xml @@ -0,0 +1,127 @@ + + + + + 0 + 3 + 1 + 0.01 + 0.50 + STOP_SIM + + 1 + result + 1 + 1 + + 100 + 0 + + 1 + 0 + 0 + + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/X0.vtp + + + + mesh/mesh-surfaces/X1.vtp + + + + mesh/mesh-surfaces/Y0.vtp + + + + mesh/mesh-surfaces/Y1.vtp + + + + mesh/mesh-surfaces/Z0.vtp + + + + mesh/mesh-surfaces/Z1.vtp + + + (0.0, 0.0, 1.0) + (1.0, 0.0, 0.0) + + + + + + + + + + true + 1 + 12 + 1e-9 + + 1.0e-2 + 1.0e5 + 0.48333 + + + 50.0 + + + + 100.0 + 8.0 + 6.0 + 12.0 + + + ST91 + + + true + true + true + true + true + true + true + true + true + true + + + + + fsils + + 1e-10 + 400 + + + + Dir + 0.0 + true + false + + + + Neu + Unsteady + load.dat + false + + + + + + + diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement.dat b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement.dat new file mode 100644 index 00000000..3ce36d1e --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:7cdbab100c4083030d6b985dee8de2aea92f8d06947f0c0c28bac88583324383 +size 1204 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement.png b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement.png new file mode 100644 index 00000000..c48ff077 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:98693767369223dd9499a8193ed2abaaadd3cf4a4d2639266b2ed2b0d0361047 +size 45457 diff --git a/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement_no_viscosity.dat b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement_no_viscosity.dat new file mode 100644 index 00000000..c6f60849 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Newtonian_viscosity/z1_face_displacement_no_viscosity.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:57946cded688a39112963bea3113a6b62e1f02de972fb4b1b7e988226e9ada87 +size 1217 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/README.md b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/README.md new file mode 100644 index 00000000..fd87ceda --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/README.md @@ -0,0 +1,27 @@ +This test case simulates pulling and releasing a slab of material +described by the Guccione material model and a pseudo-potential solid viscosity model. +In this viscosity model, the so-called viscous pseudo-potential $\Psi_{vis}$ is +given by +$$ +\Psi_{vis} = \frac{\mu}{2} \text{tr}(\dot{\mathbf{E}}) +$$ +The viscous part of the 2nd Piola-Kirchhoff stress is then given by +$$ +S_{vis} = \frac{\partial \Psi_{vis}}{\partial \dot{\mathbf{E}}} = \mu \dot{\mathbf{E}} +$$ + + + +The load profile is a 0.25s ramp to the max value, then held for 0.25s, then a 0.25s ramp down to 0, then held for 0.25s. +The load is applied on the Z1 face in the z-direction. + +![Load Profile](load.png) + +The load data is defined in `load.dat`, which can be generated with +`load.py` + +The resulting deformation is shown in the video below: +![Deformation](movie.gif) + +The z-displacement of the Z1 face is plotted below: +![Z1 Displacement](z1_face_displacement.png) diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/generate_load.py b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/generate_load.py new file mode 100644 index 00000000..88e2ac57 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/generate_load.py @@ -0,0 +1,34 @@ +import numpy as np +import os + +# Go to directory of this script +os.chdir(os.path.dirname(os.path.realpath(__file__))) + +# Number of timesteps and number of Fourier modes +n_timesteps = 1001 +n_modes = 256 + +# Generate time values from 0 to 1 +time = np.linspace(0, 1, n_timesteps) + +# Generate 0.25s ramp to max load, then hold for 0.25s, then 0.25s ramp down to 0, then hold for 0.25s +load = np.zeros(n_timesteps) +max_load= -1e2 +load[time < 0.25] = max_load * time[time < 0.25] / 0.25 +load[(time >= 0.25) & (time < 0.5)] = max_load +load[(time >= 0.5) & (time < 0.75)] = max_load * (1 - (time[(time >= 0.5) & (time < 0.75)] - 0.5) / 0.25) +load[time >= 0.75] = 0 + +# Write the time and load values to a text file +with open("load.dat", "w") as file: + file.write(f"{n_timesteps} {n_modes}\n") + for t, s in zip(time, load): + file.write(f"{t:.3f} {s:.3f}\n") + +# Plot the load values +import matplotlib.pyplot as plt +plt.plot(time, load) +plt.xlabel("Time (s)") +plt.ylabel("load (dynes/cm^2)") +plt.title("Active load") +plt.savefig("load.png") \ No newline at end of file diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/load.dat b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/load.dat new file mode 100644 index 00000000..0e638709 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/load.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1baa9ed6cf748d42b962453326733347bb3f82c1bd9c49a653568896e15e1e6c +size 13723 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/load.png b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/load.png new file mode 100644 index 00000000..5b70832a --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/load.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c58557b9830a4b793a2bb89cf130bdc4540104d8190db0c64dc59e0da1cf8d9b +size 22901 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-complete.mesh.vtu b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-complete.mesh.vtu new file mode 100644 index 00000000..f09242eb --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d74418f571657908b5c4e73f7cc5ce040fd300ac444b1c893ab23001bb6e09ae +size 7499 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/X0.vtp b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/X0.vtp new file mode 100644 index 00000000..904401ad --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/X0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3f71480fecee6f177719ab68b14bb2b61bbbce8fc745c87576fd51f3fb0b3ef8 +size 3655 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/X1.vtp b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/X1.vtp new file mode 100644 index 00000000..20b2d0e5 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/X1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:bd46bcf90705c0e062aa658846aacb9091334a56fffac81715631628cd1f097f +size 3669 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Y0.vtp b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Y0.vtp new file mode 100644 index 00000000..19c7f008 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Y0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1ef44a9a441563692f58baccf576b3c10e6df4277046d1054db5ab821d498bf0 +size 4795 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Y1.vtp b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Y1.vtp new file mode 100644 index 00000000..7ac40cb0 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Y1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4d2ad3c14f29583d9d2015b0c0ba681b66c0d2380e4b8851ca30a2cb81773183 +size 4629 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Z0.vtp b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Z0.vtp new file mode 100644 index 00000000..49112c41 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Z0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6ee78716770c3ed486e5e81d32f90226e215bc095b5cda4f5c708218b84ee407 +size 3188 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Z1.vtp b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Z1.vtp new file mode 100644 index 00000000..22122b6f --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/mesh/mesh-surfaces/Z1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c2bdc7813853e641c5a71c20f1cea215c9c0215523f51272c5650077614dcb29 +size 3245 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/movie.avi b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/movie.avi new file mode 100644 index 00000000..3a9190a5 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/movie.avi @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:a0de1f8dcee62c6bc64c602dca56942238273a19a46ac1cb78cb59e2c6c11944 +size 6604474 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/movie.gif b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/movie.gif new file mode 100644 index 00000000..e040c47c --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/movie.gif @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0d2dd0b4d57583df9fc1cd0f3b162d642e0c76cf6bb576dd5a5d9e1a2884a4ee +size 899618 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/plot_z1_face_displacement.py b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/plot_z1_face_displacement.py new file mode 100644 index 00000000..609cb008 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/plot_z1_face_displacement.py @@ -0,0 +1,145 @@ +''' +Plots the displacement of the z1 face of the mesh over time +''' + +import matplotlib.pyplot as plt +import numpy as np +import pyvista as pv +import os +import glob +import xml.etree.ElementTree as ET + +def get_timestep(result_file): + ''' + Extracts the timestep of a results_###.vtu file + ''' + # If file path is provided, get the file name + file_name = os.path.basename(result_file) + + # Get the ###.vtu part + s = file_name.split('_')[-1] + + # Get the ### part + s = s.split('.')[0] + + # Return as integer + return int(s) + +def get_start_end_step(results_folder): + """ + Automatically determine the start timestep, end timestep, and step size based on all + svFSI results file in results_folder. + + Args: + results_folder: A string of absolute file path of folder containing results of + svFSIplus simulation. This usually ends with 16-procs/ or other + number. + + Returns: + (start_timestep, end_timestep, step): A tuple of 3 integers, giving the first + timestep of results to process, the last timestep, and the step size. + + """ + + # Get list of all .vtu files in results_folder sorted by time step + list_of_files = sorted( filter( os.path.isfile, + glob.glob(os.path.join(results_folder, '*.vtu')) ), key = get_timestep) + + # Get start time from the first result file (list_of_files[0]) + start_file_name = os.path.basename(list_of_files[0]) + start_timestep = int("".join([i for i in start_file_name if i.isdigit()])) + print('Start timestep:', start_timestep) + + # Get end time from the last result file (list_of_files[-1]) + end_file_name = os.path.basename(list_of_files[-1]) + end_timestep = int("".join([i for i in end_file_name if i.isdigit()])) + print('End timestep:', end_timestep) + + # Get step size by looking at second time step + start_plus_one_file_name = os.path.basename(list_of_files[1]) + start_time_plus_one = int("".join([i for i in start_plus_one_file_name if i.isdigit()])) + step = start_time_plus_one - start_timestep + print('Step:', step) + + return (start_timestep, end_timestep, step) + +def get_timestep_size(svfsiplus_xml_file): + ''' + Reads the timestep size from an svFSIplus XML file. + + Args: + svfsiplus_xml_file: The file path to the svFSIplus XML input file + + Returns: + timestep_size: The size of the timestep in the input file + ''' + # Parse the XML file + tree = ET.parse(svfsiplus_xml_file) + root = tree.getroot() + + # Find the element containing the Time_step_size + timestep_size = None + for elem in root.iter(): + if elem.tag == 'Time_step_size': + timestep_size = float(elem.text) + break + + return timestep_size + + +# Move to this script's directory +os.chdir(os.path.dirname(os.path.realpath(__file__))) + +# Load the z1 face +mesh = pv.read("mesh/mesh-surfaces/Z1.vtp") + +# Get the start timestep, end timestep, and step size +start_timestep, end_timestep, step = get_start_end_step("4-procs/") + +# Get the timestep size +timestep_size = get_timestep_size("svFSIPlus.xml") + +# Get the time values +time = np.arange(start_timestep, end_timestep + 1, step) * timestep_size + +# Loop over the results files +displacements = [] +for n in range(start_timestep, end_timestep + 1, step): + # Load the results file + result_file = f"4-procs/result_{n:03d}.vtu" + result = pv.read(result_file) + + # Sample displacement at the z1 face + displacement = mesh.sample(result).point_data['Displacement'][:, 2].mean() + + # Append to the list + displacements.append(displacement) + +# Save the displacements to a text file +with open("z1_face_displacement.dat", "w") as file: + file.write(f"{len(time)}\n") + for t, s in zip(time, displacements): + file.write(f"{t:.3f} {s:.3f}\n") + + +# Plot the displacement +plt.plot(time, displacements, label="With viscosity") +plt.xlabel("Time (s)") +plt.ylabel("z Displacement (cm)") +plt.title("Z1 face displacement") + +# Plot the displacement with no viscosity +if os.path.exists("z1_face_displacement_no_viscosity.dat"): + with open("z1_face_displacement_no_viscosity.dat", "r") as file: + lines = file.readlines() + n = int(lines[0]) + time_no_viscosity = np.zeros(n) + displacements_no_viscosity = np.zeros(n) + for i, line in enumerate(lines[1:]): + time_no_viscosity[i], displacements_no_viscosity[i] = map(float, line.split()) + + plt.plot(time_no_viscosity, displacements_no_viscosity, label="Without viscosity") + + +plt.legend() +plt.savefig("z1_face_displacement.png") \ No newline at end of file diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/result_001.vtu b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/result_001.vtu new file mode 100644 index 00000000..d0a09cf3 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/result_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:dafbffac82bd8ed8f5945656a75ee05bf61c73b7adaed04c49c730bf75223687 +size 81248 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/svFSIplus.xml b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/svFSIplus.xml new file mode 100644 index 00000000..28a7309d --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/svFSIplus.xml @@ -0,0 +1,127 @@ + + + + + 0 + 3 + 1 + 0.01 + 0.50 + STOP_SIM + + 1 + result + 1 + 1 + + 100 + 0 + + 1 + 0 + 0 + + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/X0.vtp + + + + mesh/mesh-surfaces/X1.vtp + + + + mesh/mesh-surfaces/Y0.vtp + + + + mesh/mesh-surfaces/Y1.vtp + + + + mesh/mesh-surfaces/Z0.vtp + + + + mesh/mesh-surfaces/Z1.vtp + + + (0.0, 0.0, 1.0) + (1.0, 0.0, 0.0) + + + + + + + + + + true + 1 + 12 + 1e-9 + + 1.0e-2 + 1.0e5 + 0.48333 + + + 50.0 + + + + 100.0 + 8.0 + 6.0 + 12.0 + + + ST91 + + + true + true + true + true + true + true + true + true + true + true + + + + + fsils + + 1e-10 + 400 + + + + Dir + 0.0 + true + false + + + + Neu + Unsteady + load.dat + false + + + + + + + diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/z1_face_displacement.dat b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/z1_face_displacement.dat new file mode 100644 index 00000000..52145635 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/z1_face_displacement.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:900e23f9b4d800bf654699775bf1d91ce31b4d7e5807c79702bbd65de967463a +size 1204 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/z1_face_displacement.png b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/z1_face_displacement.png new file mode 100644 index 00000000..ccd3b05b --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/z1_face_displacement.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ab9c4afaa95073f614ee06d3523caa606ad422d79b0e1a6c61567c4b975f289e +size 44709 diff --git a/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/z1_face_displacement_no_viscosity.dat b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/z1_face_displacement_no_viscosity.dat new file mode 100644 index 00000000..c6f60849 --- /dev/null +++ b/tests/cases/ustruct/tensile_adventitia_Potential_viscosity/z1_face_displacement_no_viscosity.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:57946cded688a39112963bea3113a6b62e1f02de972fb4b1b7e988226e9ada87 +size 1217 diff --git a/tests/test_struct.py b/tests/test_struct.py index 4fda6376..9b02081a 100644 --- a/tests/test_struct.py +++ b/tests/test_struct.py @@ -68,3 +68,11 @@ def test_LV_NeoHookean_passive_sv0D(n_proc): def test_tensile_adventitia_Guccione_active(n_proc): test_folder = "tensile_adventitia_Guccione_active" run_with_reference(base_folder, test_folder, fields, n_proc, t_max=2) + +def test_tensile_adventitia_Newtonian_viscosity(n_proc): + test_folder = "tensile_adventitia_Newtonian_viscosity" + run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) + +def test_tensile_adventitia_Potential_viscosity(n_proc): + test_folder = "tensile_adventitia_Potential_viscosity" + run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) diff --git a/tests/test_ustruct.py b/tests/test_ustruct.py index f56a5fab..f3fb9112 100644 --- a/tests/test_ustruct.py +++ b/tests/test_ustruct.py @@ -57,3 +57,11 @@ def test_LV_NeoHookean_passive_sv0D(n_proc): test_folder = "LV_NeoHookean_passive_sv0D" run_with_reference(base_folder, test_folder, fields, n_proc, t_max=3) + +def test_tensile_adventitia_Newtonian_viscosity(n_proc): + test_folder = "tensile_adventitia_Newtonian_viscosity" + run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) + +def test_tensile_adventitia_Potential_viscosity(n_proc): + test_folder = "tensile_adventitia_Potential_viscosity" + run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) \ No newline at end of file