Check constitutive model ustruct 320 (SimVascular#321)

* Change file name.

* Add functions to check for a valid constitutive model for ustruct equations.

Author: ktbolt

Code/Source/solver/Parameters.cpp

@@ -63,6 +63,7 @@
 #include "Parameters.h"
 #include "consts.h"
 #include "LinearAlgebra.h"
+#include "ustruct.h"
 
 #include <iostream>
 #include <regex>
@@ -791,6 +792,20 @@ void ConstitutiveModelParameters::set_values(tinyxml2::XMLElement* xml_elem)
   value_set = true;
 }
 
+/// @brief Check if a constitutive model is valid for the given equation.
+//
+void ConstitutiveModelParameters::check_constitutive_model(const Parameter<std::string>& eq_type_str)
+{
+  auto eq_type = consts::equation_name_to_type.at(eq_type_str.value());
+  auto model = consts::constitutive_model_name_to_type.at(type.value());
+
+  if (eq_type == consts::EquationType::phys_ustruct) {
+    if (! ustruct::constitutive_model_is_valid(model)) {
+      throw std::runtime_error("The " + type.value() + " constitutive model is not valid for ustruct equations.");
+    }
+  }
+}
+
 //////////////////////////////////////////////////////////
 //                  CoupleCplBCParameters               //
 //////////////////////////////////////////////////////////
@@ -1864,7 +1879,6 @@ void EquationParameters::set_values(tinyxml2::XMLElement* eq_elem)
   //
   while (item != nullptr) {
     auto name = std::string(item->Value());
-    //std::cout << "[EquationParameters::set_values] name: " << name << std::endl;
 
     if (name == BodyForceParameters::xml_element_name_) {
       auto bf_params = new BodyForceParameters();
@@ -1878,6 +1892,7 @@ void EquationParameters::set_values(tinyxml2::XMLElement* eq_elem)
 
     } else if (name == ConstitutiveModelParameters::xml_element_name_) {
       default_domain->constitutive_model.set_values(item);
+      default_domain->constitutive_model.check_constitutive_model(type);
 
     } else if (name == CoupleCplBCParameters::xml_element_name_) {
       couple_to_cplBC.set_values(item);

Code/Source/solver/Parameters.h

@@ -510,6 +510,7 @@ class ConstitutiveModelParameters : public ParameterLists
   public:
     ConstitutiveModelParameters();
     void print_parameters();
+    void check_constitutive_model(const Parameter<std::string>& eq_type);
     bool defined() const { return value_set; };
     void set_values(tinyxml2::XMLElement* modl_params);
     static const std::string xml_element_name_;

Code/Source/solver/mat_models.cpp

@@ -408,21 +408,13 @@ void compute_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType&
   // Now, add isochoric and total stress, elasticity tensors
   switch (stM.isoType) {
     case ConstitutiveModelType::stIso_lin: {
-      if (ustruct) {
-        throw std::runtime_error("[compute_pk2cc] Linear isotropic material model not valid for ustruct physics.");
-      }
-
       double g1 = stM.C10;    // mu
       S += g1*Idm;
       return; 
     } break;
 
     // St.Venant-Kirchhoff
     case ConstitutiveModelType::stIso_StVK: {
-      if (ustruct) {
-        throw std::runtime_error("[compute_pk2cc] St.Venant-Kirchhoff material model not valid for ustruct physics.");
-      }
-      
       double g1 = stM.C10;         // lambda
       double g2 = stM.C01 * 2.0;   // 2*mu
 
@@ -432,10 +424,6 @@ void compute_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType&
 
     // modified St.Venant-Kirchhoff
     case ConstitutiveModelType::stIso_mStVK: {
-      if (ustruct) {
-        throw std::runtime_error("[compute_pk2cc] Modified St.Venant-Kirchhoff material model not valid for ustruct physics.");
-      }
-
       double g1 = stM.C10; // kappa
       double g2 = stM.C01;  // mu
 

Code/Source/solver/ustruct.cpp

@@ -210,6 +210,21 @@ void b_ustruct_3d(const ComMod& com_mod, const int eNoN, const double w, const V
   }
 }
 
+/// @brief Check is a constitutive model is valid for the ustruct equation.
+//
+bool constitutive_model_is_valid(consts::ConstitutiveModelType model)
+{ 
+  using namespace consts;
+    
+  static std::set<ConstitutiveModelType> unsupported_models { 
+    ConstitutiveModelType::stIso_lin,
+    ConstitutiveModelType::stIso_StVK,
+    ConstitutiveModelType::stIso_mStVK
+  };
+
+  return unsupported_models.count(model) == 0;
+}   
+
 /// @brief Reproduces Fortran CONSTRUCT_uSOLID.
 //
 void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array<double>& Ag, 

Code/Source/solver/ustruct.h

@@ -43,6 +43,8 @@ void b_ustruct_3d(const ComMod& com_mod, const int eNoN, const double w, const V
     const Array<double>& Nx, const Array<double>& dl, const Vector<double>& hl, const Vector<double>& nV, 
     Array<double>& lR, Array3<double>& lK, Array3<double>& lKd);
 
+bool constitutive_model_is_valid(consts::ConstitutiveModelType model);
+
 void construct_usolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array<double>& Ag, const Array<double>& Yg, 
     const Array<double>& Dg);