diff --git a/Code/Source/README.md b/Code/Source/README.md
index bbe377b7..f1c5bb17 100644
--- a/Code/Source/README.md
+++ b/Code/Source/README.md
@@ -13,6 +13,9 @@ This document describes some of the implementation details of svFSIplus C++ code
[Performance and Accuracy](#performance)
[Potential Problems](#problems)
[Implementation Details](#cpp_programming)
+[Coding Standards](#coding_standards)
+[Coding Guidelines](#coding_guidelines)
+
@@ -790,6 +793,341 @@ which are used like this
if ((eq.phys == Equation_CMM) && com_mod.cmmInit) {
```
+
+
+
+
+
Coding Standards
+
+This section describes the coding standards and guidelines that must be followed when adding new code to svFSIplus.
+
+Coding standards are important for producing good software
+
+- Gives a uniform appearance to the code written by different developers
+- Improves readability and maintainability
+- Helps in code reuse and to detect errors
+- Promotes sound programming practices and increases programmer efficiency
+
+Code that does not follow coding standards will be rejected during code review. However, program structures not mentioned in the
+following sections may be coded in any (non-hideous) manner favoured by the developer.
+
+Note that svFSIplus maintained the program structure and naming conventions used by the svFSI Fortran code so some of the
+following coding standards may be violated.
+
+ Indentation
+
+Indentation refers to the spaces at the beginning of a code line. It is a fundamental aspect of code styling and improves
+readability by showing the overall structure of the code.
+
+Indentation is two spaces for all programming structures: functions, loops, if-then blocks, etc. Do not use tabs to indent.
+
+The if-else class of statements should have the following form
+```
+if (condition) {
+ ...
+}
+
+if (condition) {
+ ...
+} else {
+ ...
+}
+```
+
+A for statement should have the following form
+```
+for (initialization; condition; update) {
+ statements;
+}
+```
+
+A switch statement should have the following form
+```
+switch (condition) {
+ case ABC:
+ statements;
+ break;
+
+ default:
+ statements;
+ break;
+}
+```
+
+The braces indicating a function body are placed in column 1, function body statements are indented by two spaces.
+```
+void find_max(const int n, const double values[])
+{
+ int num_values{0};
+
+ for (int i = 0; i < n; i++) {
+ if (values[i] != 0.0) {
+ num_values += 1;
+ }
+ }
+}
+```
+
+ White Space
+
+Whitespace is a term that refers to the spaces and newlines that are used to make code more readable.
+
+The following white space conventions should be followed
+- Assignment operators always have spaces around them.
+- Other binary operators usually have spaces around them, but it’s OK to remove spaces around *, / and %.
+- Parentheses should have no internal padding.
+- No spaces separating unary operators and their arguments.
+- C++ reserved words should be followed by a white space.
+- Commas should be followed by a white space.
+- Colons should be surrounded by white space.
+- Semicolons in for statements should be followed by a space character.
+
+```
+v = w * x + y / z; // good
+v = w*x + y/z; // also ok
+v = w * (x + z); // good
+while (true) { // bad while(true) ...
+doSomething(a, b, c, d); // bad doSomething(a,b,c,d);
+for (i = 0; i < 10; i++) { // bad for (i=0;i<10;i++){
+```
+
+Some complex expressions may be better organized without single separating spaces. The following could be written using
+spaces between sub-expressions only
+```
+double Inv2 = 0.50 * (Inv1*Inv1 - J4d*mat_trace(mat_mul(C,C), nsd));
+```
+or written using double spacing between sub-expressions
+```
+double Inv2 = 0.50 * (Inv1 * Inv1 - J4d * mat_trace(mat_mul(C,C), nsd));
+```
+
+Use newlines often to separate logical blocks of code: for-loops, if statements, related code blocks.
+```
+ for (int e = 0; e < lM.nEl; e++) {
+ int iDmn = all_fun::domain(com_mod, lM, cEq, e);
+ auto cPhys = com_mod.eq[cEq].dmn[iDmn].phys;
+
+ for (int a = 0; a < lM.eNoN; a++) {
+ int Ac = lM.IEN(a,e);
+ xl.set_col(a, com_mod.x.col(Ac));
+
+ if (com_mod.mvMsh) {
+ for (int i = 0; i < nsd; i++) {
+ dol(i,a) = com_mod.Do(i + nsd + 1, Ac);
+ }
+ }
+ }
+
+ if (com_mod.mvMsh) {
+ xl = xl + dol;
+ }
+```
+
+ Namespaces
+
+Refer to the elements in the std namespace by explicit qualification using std::.
+```
+std::cout << "Use std namespace." << std::endl;
+std::string file_name;
+std::vector outS(nOut+1);
+```
+
+It is acceptable to use unqualified names for svFSIplus namespaces
+```
+using namespace consts;
+double dmp = dmn.prop.at(PhysicalProperyType::damping);
+```
+
+ Naming Conventions
+
+Program readability is improved by using names that would be clear to any developer.
+
+- Use names that describe purpose or intent
+- There's no reason to be terse when naming things. It is more important to make code immediately
+ understandable by a new reader
+- Minimize the use of specialized abbreviations (e.g., those from a paper or text book) that may
+ in general be unfamiliar to other developers
+- Do not abbreviate names especially for class methods/data members
+ - command instead of cmd
+ - initialize instead of init
+ - compute_average() instead of comp_avg()
+- Local names may be abbreviated for local variables when context is clear
+ - num_points is OK instead of number_of_points
+- Variables with a large scope should have long names, variables with a small scope may have short names
+- Names for boolean variables and methods should be obviously boolean
+ - use is_valid instead of flag
+
+
+
+ Styles
+
+Two naming styles are used svFSIplus
+
+- CamelCase
+
+- snake_case
+
+CamelCase is a way to separate the words in a phrase by making the first letter of each word capitalized
+and not using spaces.
+
+
+ File Names
+
+C++ files should end in .cpp and header files should end in .h.
+
+Files that contain a single class should have the name of the class, including capitalization.
+- MyClass.h, MyClass.cc
+
+
+ Type Names
+
+Type names are in CamelCase.
+
+```
+class Simulation { }
+CepMod cep_mod;
+```
+
+ Variable and Functions Names
+
+Variable and function names use snake_case.
+```
+std::string history_file_name;
+void b_ustruct_3d(const ComMod& com_mod);
+```
+Data members of classes additionally have trailing underscores.
+
+
+
+Comments are absolutely vital to keeping code readable. While comments are very important, the best code is self-documenting.
+Giving sensible names to types and variables is much better than using obscure names that you must then explain through comments.
+
+Don't literally describe what code does unless the behavior is nonobvious to a reader who understands C++ well.
+Instead, provide higher level comments that describe why the code does what it does, or make the code self describing.
+
+Comments should be included relative to their position in the code
+```
+// yes
+while (true) {
+ // Do something
+ something();
+}
+
+// no
+while (true) {
+// Do something now
+ something();
+}
+```
+
+ File Comments
+
+Start each file with license boilerplate. Every file should contain license boilerplate.
+
+File comments describe the contents of a file. A .h file will contain comments describing any classes defined there. A .cpp file
+will contain comments describing the purpose of the functions or class methods defined there.
+
+Do not duplicate comments in both the .h and the .cpp files.
+
+ Class Comments
+
+Every non-obvious class or struct declaration should have an accompanying comment that describes what it is for and how it should be used.
+
+Comments describing the use of the class should go together with its interface definition; comments about the class operation and
+implementation should accompany the implementation of the class's methods.
+
+ Function Comments
+
+Every function declaration should have comments immediately preceding it that describe what the function does and how to use it. If there
+is anything tricky about how a function does its job, the function definition should have an explanatory comment.
+
+The function implementation should have comments describing tricky, non-obvious, interesting, or important parts of the code.
+
+Function comments should follow Doxygen format for API functions.
+
+Non API functions should have the form
+```
+Parameters
+ param1[in] Description
+ param1[out] Description
+
+Returns
+```
+
+
+
+ Variables
+
+Variables should be initialized where they are declared when possible. This ensures that variables are valid at any time.
+```
+double average{0.0};
+int num_points = 0;
+```
+
+Variables must never have dual meaning. This ensures that all concepts are represented uniquely.
+
+Global variable use should be minimized. In C++ there is no reason that global variables need to be used at all.
+
+Variables should be declared in the smallest scope possible. By keeping the operations on a variable within a small scope
+it is easier to control the effects and side effects of the variable.
+
+
+ General Programming
+
+Use nullptr instead of 0 and NULL.
+
+Use const rather than #define statements.
+```
+// Replace
+#define A_POWER_OF_TWO 16
+
+// with
+int const A_POWER_OF_TWO = 16;
+```
+
+Avoid deeply nested code. Code that is too deeply nested is hard to both read and debug.
+One should replace excessive nesting with function calls.
+
+
+
+
+
+ Coding Guidelines
+
+This section describes the coding guidelines that are recommend when adding new code to svFSIplus.
+
+ Enums
+
+Where possible, put enums in appropriate classes
+```
+class Grade {
+ enum { HIGH, MIDDLE, LOW };
+
+ Grade() {}
+ ...
+};
+```
+
+ Type Conversions
+
+Type conversions should be avoided if possible.
+
+When required, type conversions must always be done explicitly using C++ style casts. Never rely on implicit type conversion.
+
+```
+double r = static_cast(i) / 3.0;
+```
+
+ Function Parameters
+
+Arguments that are non-primitive types and will not be modified should be passed by const reference.
+```
+void calc_elem_ar(const ComMod& com_mod, const CmMod& cm_mod, mshType& lM, bool& rflag)
+```
+
+Output parameters should grouped at the end of the function's parameters.
+
+
diff --git a/Code/Source/svFSI/Array.h b/Code/Source/svFSI/Array.h
index e3145a22..ed562591 100644
--- a/Code/Source/svFSI/Array.h
+++ b/Code/Source/svFSI/Array.h
@@ -907,6 +907,21 @@ class Array
return *this;
}
+ //--------
+ // negate
+ //--------
+ //
+ Array operator-() const
+ {
+ Array result(nrows_, ncols_);
+ for (int j = 0; j < ncols_; j++) {
+ for (int i = 0; i < nrows_; i++) {
+ result(i,j) = -(data_[i + j*nrows_]);
+ }
+ }
+ return result;
+ }
+
// Compound subtract assignment.
//
Array operator-=(const T value) const
@@ -996,7 +1011,7 @@ class Array
return sum;
}
- T* data() {
+ T* data() const {
return data_;
}
diff --git a/Code/Source/svFSI/Array3.h b/Code/Source/svFSI/Array3.h
index 8619c262..67816dbe 100644
--- a/Code/Source/svFSI/Array3.h
+++ b/Code/Source/svFSI/Array3.h
@@ -296,6 +296,30 @@ class Array3
allocate(num_rows, num_cols, num_slices);
}
+ //------------
+ // set_values
+ //------------
+ // Set the array values from an Array with the equivalent
+ // number of values.
+ //
+ // This sort of replicates the Fortan reshape function.
+ //
+ void set_values(Array& rhs)
+ {
+ int rhs_size = rhs.size();
+
+ if (size_ != rhs_size) {
+ throw std::runtime_error("The RHS size " + std::to_string(rhs_size) + " does not have the same size " +
+ std::to_string(size_) + " of this array.");
+ }
+
+ auto rhs_data = rhs.data();
+
+ for (int i = 0; i < size_; i++) {
+ data_[i] = rhs_data[i];
+ }
+ }
+
//------
// read
//------
diff --git a/Code/Source/svFSI/ComMod.h b/Code/Source/svFSI/ComMod.h
index ee2981e2..03f18c91 100644
--- a/Code/Source/svFSI/ComMod.h
+++ b/Code/Source/svFSI/ComMod.h
@@ -364,6 +364,15 @@ class stModelType
// Collagen fiber dispersion parameter (Holzapfel model)
double kap = 0.0;
+ // Heaviside function parameter (Holzapfel-Ogden model)
+ double khs = 100.0;
+
+ // Lee-Sacks model
+ double a0 = 0.0;
+ double b1 = 0.0;
+ double b2 = 0.0;
+ double mu0 = 0.0;
+
// Fiber reinforcement stress
fibStrsType Tf;
};
@@ -666,7 +675,7 @@ class cntctModelType
{
public:
// Contact model
- int cType = 0;
+ consts::ContactModelType cType = consts::ContactModelType::cntctM_NA;
// Penalty parameter
double k = 0.0;
diff --git a/Code/Source/svFSI/Parameters.cpp b/Code/Source/svFSI/Parameters.cpp
index 8dbb5674..78454a9b 100644
--- a/Code/Source/svFSI/Parameters.cpp
+++ b/Code/Source/svFSI/Parameters.cpp
@@ -139,6 +139,9 @@ void Parameters::read_xml(std::string file_name)
// Get general parameters.
general_simulation_parameters.set_values(root_element);
+ // Set Contact values.
+ set_contact_values(root_element);
+
// Set Add_mesh values.
set_mesh_values(root_element);
@@ -149,6 +152,21 @@ void Parameters::read_xml(std::string file_name)
set_equation_values(root_element);
}
+//--------------------
+// set_contact_values
+//--------------------
+//
+void Parameters::set_contact_values(tinyxml2::XMLElement* root_element)
+{
+ auto item = root_element->FirstChildElement(ContactParameters::xml_element_name_.c_str());
+
+ if (item == nullptr) {
+ return;
+ }
+
+ contact_parameters.set_values(item);
+}
+
//---------------------
// set_equation_values
//---------------------
@@ -385,7 +403,7 @@ BoundaryConditionParameters::BoundaryConditionParameters()
set_parameter("Ramp_function", false, !required, ramp_function);
- set_parameter("Shell_bc_type", "", !required, shell_bc_type);
+ set_parameter("CST_shell_bc_type", "", !required, cst_shell_bc_type);
set_parameter("Spatial_profile_file_path", "", !required, spatial_profile_file_path);
set_parameter("Spatial_values_file_path", "", !required, spatial_values_file_path);
set_parameter("Stiffness", 1.0, !required, stiffness);
@@ -473,6 +491,7 @@ const std::string ConstitutiveModelParameters::xml_element_name_ = "Constitutive
// [TODO] Should use the types defined in consts.h.
const std::string ConstitutiveModelParameters::GUCCIONE_MODEL = "Guccione";
const std::string ConstitutiveModelParameters::HGO_MODEL = "HGO";
+const std::string ConstitutiveModelParameters::LEE_SACKS = "Lee-Sacks";
const std::string ConstitutiveModelParameters::NEOHOOKEAN_MODEL = "neoHookean";
const std::string ConstitutiveModelParameters::STVENANT_KIRCHHOFF_MODEL = "stVenantKirchhoff";
@@ -487,6 +506,8 @@ const std::map ConstitutiveModelParameters::constituti
{ConstitutiveModelParameters::HGO_MODEL, ConstitutiveModelParameters::HGO_MODEL},
+ {ConstitutiveModelParameters::LEE_SACKS, ConstitutiveModelParameters::LEE_SACKS},
+
{ConstitutiveModelParameters::NEOHOOKEAN_MODEL, ConstitutiveModelParameters::NEOHOOKEAN_MODEL},
{"nHK", ConstitutiveModelParameters::NEOHOOKEAN_MODEL},
@@ -503,6 +524,7 @@ using SetConstitutiveModelParamMapType = std::map void {cp->guccione.set_values(params);}},
{ConstitutiveModelParameters::HGO_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->holzapfel_gasser_ogden.set_values(params);}},
+ {ConstitutiveModelParameters::LEE_SACKS, [](CmpType cp, CmpXmlType params) -> void {cp->lee_sacks.set_values(params);}},
{ConstitutiveModelParameters::NEOHOOKEAN_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->neo_hookean.set_values(params);}},
{ConstitutiveModelParameters::STVENANT_KIRCHHOFF_MODEL, [](CmpType cp, CmpXmlType params) -> void {cp->stvenant_kirchhoff.set_values(params);}},
};
@@ -514,10 +536,12 @@ using PrintConstitutiveModelParamMapType = std::map void {cp->guccione.print_parameters();}},
{ConstitutiveModelParameters::HGO_MODEL, [](CmpType cp) -> void {cp->holzapfel_gasser_ogden.print_parameters();}},
+ {ConstitutiveModelParameters::LEE_SACKS, [](CmpType cp) -> void {cp->lee_sacks.print_parameters();}},
{ConstitutiveModelParameters::NEOHOOKEAN_MODEL, [](CmpType cp) -> void {cp->neo_hookean.print_parameters();}},
{ConstitutiveModelParameters::STVENANT_KIRCHHOFF_MODEL, [](CmpType cp) -> void {cp->stvenant_kirchhoff.print_parameters();}},
};
+
//-------------------------------------
// ConstitutiveModelGuccioneParameters
//-------------------------------------
@@ -642,6 +666,47 @@ void HolzapfelGasserOgdenParameters::print_parameters()
}
}
+//--------------------
+// LeeSacksParameters
+//--------------------
+//
+LeeSacksParameters::LeeSacksParameters()
+{
+ // A parameter that must be defined.
+ bool required = true;
+
+ set_parameter("a", 0.0, required, a);
+ set_parameter("a0", 0.0, required, a);
+ set_parameter("b1", 0.0, required, b1);
+ set_parameter("b2", 0.0, required, b2);
+ set_parameter("mu0", 0.0, required, mu0);
+
+ set_xml_element_name("Constitutive_model type=Lee-Sacks");
+}
+
+void LeeSacksParameters::set_values(tinyxml2::XMLElement* xml_elem)
+{
+ std::string error_msg = "Unknown Constitutive_model type=Lee-Sacks XML element '";
+
+ using std::placeholders::_1;
+ using std::placeholders::_2;
+ std::function ftpr =
+ std::bind( &LeeSacksParameters::set_parameter_value, *this, _1, _2);
+
+ xml_util_set_parameters(ftpr, xml_elem, error_msg);
+
+ value_set = true;
+}
+
+void LeeSacksParameters::print_parameters()
+{
+ std::cout << "Lee-Sacks: " << std::endl;
+ auto params_name_value = get_parameter_list();
+ for (auto& [ key, value ] : params_name_value) {
+ std::cout << key << ": " << value << std::endl;
+ }
+}
+
//------------------------
// MooneyRivlinParameters
//------------------------
@@ -1304,7 +1369,7 @@ DomainParameters::DomainParameters()
set_parameter("Poisson_ratio", 0.3, !required, poisson_ratio);
set_parameter("Relative_tolerance", 1e-4, !required, relative_tolerance);
- set_parameter("Shell_thicknes", 0.0, !required, shell_thickness);
+ set_parameter("Shell_thickness", 0.0, !required, shell_thickness);
set_parameter("Solid_density", 0.5, !required, solid_density);
set_parameter("Solid_viscosity", 0.9, !required, solid_viscosity);
set_parameter("Source_term", 0.0, !required, source_term);
@@ -1599,6 +1664,77 @@ void ECGLeadsParameters::print_parameters()
}
}
+//////////////////////////////////////////////////////////
+// ContactParameters //
+//////////////////////////////////////////////////////////
+
+// Process parameters for the 'Contact' XML element
+// used to specify parameters for contact computation.
+
+// Define the XML element name for contact parameters.
+const std::string ContactParameters::xml_element_name_ = "Contact";
+
+//--------------------
+// EquationParameters
+//--------------------
+//
+ContactParameters::ContactParameters()
+{
+ set_xml_element_name(xml_element_name_);
+
+ // A parameter that must be defined.
+ bool required = true;
+
+ // Contact model.
+ model = Parameter("model", "", required);
+
+ // Define contact parameters.
+ //
+ set_parameter("Closest_gap_to_activate_penalty", 1.0, !required, closest_gap_to_activate_penalty);
+ set_parameter("Desired_separation", 0.05, !required, desired_separation);
+ set_parameter("Min_norm_of_face_normals", 0.7, !required, min_norm_of_face_normals);
+ set_parameter("Penalty_constant", 1e5, !required, penalty_constant);
+}
+
+void ContactParameters::print_parameters()
+{
+ std::cout << std::endl;
+ std::cout << "-------------------" << std::endl;
+ std::cout << "Contact Parameters" << std::endl;
+ std::cout << "-------------------" << std::endl;
+ std::cout << model.name() << ": " << model.value() << std::endl;
+
+ auto params_name_value = get_parameter_list();
+ for (auto& [ key, value ] : params_name_value) {
+ std::cout << key << ": " << value << std::endl;
+ }
+}
+
+//------------
+// set_values
+//------------
+void ContactParameters::set_values(tinyxml2::XMLElement* xml_elem)
+{
+ using namespace tinyxml2;
+ std::string error_msg = "Unknown " + xml_element_name_ + " XML element '";
+
+ // Get the 'type' from the element.
+ const char* mname;
+ auto result = xml_elem->QueryStringAttribute("model", &mname);
+ if (mname == nullptr) {
+ throw std::runtime_error("No MODEL given in the XML element.");
+ }
+ model.set(std::string(mname));
+
+ using std::placeholders::_1;
+ using std::placeholders::_2;
+
+ std::function ftpr =
+ std::bind( &ProjectionParameters::set_parameter_value, *this, _1, _2);
+
+ xml_util_set_parameters(ftpr, xml_elem, error_msg);
+}
+
//////////////////////////////////////////////////////////
// EquationParameters //
//////////////////////////////////////////////////////////
diff --git a/Code/Source/svFSI/Parameters.h b/Code/Source/svFSI/Parameters.h
index d5caf8e8..737d9ab8 100644
--- a/Code/Source/svFSI/Parameters.h
+++ b/Code/Source/svFSI/Parameters.h
@@ -439,6 +439,25 @@ class ParameterLists
// The following classes are used to store parameters for
// various constitutive models.
+//--------------------
+// LeeSacksParameters
+//--------------------
+//
+class LeeSacksParameters : public ParameterLists
+{
+ public:
+ LeeSacksParameters();
+ bool defined() const { return value_set; };
+ void set_values(tinyxml2::XMLElement* con_model_params);
+ void print_parameters();
+ Parameter a;
+ Parameter a0;
+ Parameter b1;
+ Parameter b2;
+ Parameter mu0;
+ bool value_set = false;
+};
+
//--------------------
// GuccioneParameters
//--------------------
@@ -548,6 +567,7 @@ class ConstitutiveModelParameters : public ParameterLists
// Model types supported.
static const std::string GUCCIONE_MODEL;
static const std::string HGO_MODEL;
+ static const std::string LEE_SACKS;
static const std::string NEOHOOKEAN_MODEL;
static const std::string STVENANT_KIRCHHOFF_MODEL;
static const std::map constitutive_model_types;
@@ -558,6 +578,7 @@ class ConstitutiveModelParameters : public ParameterLists
GuccioneParameters guccione;
HolzapfelParameters holzapfel;
HolzapfelGasserOgdenParameters holzapfel_gasser_ogden;
+ LeeSacksParameters lee_sacks;
MooneyRivlinParameters mooney_rivlin;
NeoHookeanParameters neo_hookean;
StVenantKirchhoffParameters stvenant_kirchhoff;
@@ -716,7 +737,7 @@ class BoundaryConditionParameters : public ParameterLists
Parameter profile;
Parameter ramp_function;
- Parameter shell_bc_type;
+ Parameter cst_shell_bc_type;
Parameter spatial_profile_file_path;
Parameter spatial_values_file_path;
Parameter stiffness;
@@ -1086,6 +1107,34 @@ class RemesherParameters : public ParameterLists
Parameter frequency_for_copying_data;
};
+//-------------------
+// ContactParameters
+//-------------------
+// The ContactParameters class stores parameters for the 'Contact''
+// XML element used to specify parameter values for contact
+// computations.
+//
+class ContactParameters : public ParameterLists
+{
+ public:
+ ContactParameters();
+
+ static const std::string xml_element_name_;
+
+ void print_parameters();
+ void set_values(tinyxml2::XMLElement* xml_elem);
+
+ Parameter closest_gap_to_activate_penalty;
+
+ Parameter desired_separation;
+
+ Parameter min_norm_of_face_normals;
+
+ Parameter model;
+
+ Parameter penalty_constant;
+};
+
//--------------------
// EquationParameters
//--------------------
@@ -1291,11 +1340,13 @@ class Parameters {
void print_parameters();
void read_xml(std::string file_name);
+ void set_contact_values(tinyxml2::XMLElement* root_element);
void set_equation_values(tinyxml2::XMLElement* root_element);
void set_mesh_values(tinyxml2::XMLElement* root_element);
void set_projection_values(tinyxml2::XMLElement* root_element);
// Objects representing each parameter section of XML file.
+ ContactParameters contact_parameters;
GeneralSimulationParameters general_simulation_parameters;
std::vector mesh_parameters;
std::vector equation_parameters;
diff --git a/Code/Source/svFSI/VtkData.cpp b/Code/Source/svFSI/VtkData.cpp
index 3221a3a9..bf61687d 100644
--- a/Code/Source/svFSI/VtkData.cpp
+++ b/Code/Source/svFSI/VtkData.cpp
@@ -291,10 +291,13 @@ void VtkVtuData::VtkVtuDataImpl::set_connectivity(const int nsd, const ArrayGetPoints()->GetNumberOfPoints();
unsigned char vtk_cell_type;
- //std::cout << "[VtkVtuData.set_connectivity] " << std::endl;
- //std::cout << "[VtkVtuData.set_connectivity] num_elems: " << num_elems << std::endl;
- //std::cout << "[VtkVtuData.set_connectivity] np_elem: " << np_elem << std::endl;
- //std::cout << "[VtkVtuData.set_connectivity] num_coords: " << num_coords << std::endl;
+ /*
+ std::cout << "[VtkVtuData.set_connectivity] " << std::endl;
+ std::cout << "[VtkVtuData.set_connectivity] nsd: " << nsd << std::endl;
+ std::cout << "[VtkVtuData.set_connectivity] num_elems: " << num_elems << std::endl;
+ std::cout << "[VtkVtuData.set_connectivity] np_elem: " << np_elem << std::endl;
+ std::cout << "[VtkVtuData.set_connectivity] num_coords: " << num_coords << std::endl;
+ */
if (nsd == 2) {
@@ -314,7 +317,10 @@ void VtkVtuData::VtkVtuDataImpl::set_connectivity(const int nsd, const Array& Dg)
//
void set_bf_l(ComMod& com_mod, bfType& lBf, mshType& lM, const Array& Dg)
{
- #define debug_set_bf_l
+ #define n_debug_set_bf_l
auto& cm = com_mod.cm;
#ifdef debug_set_bf_l
DebugMsg dmsg(__func__, com_mod.cm.idcm());
diff --git a/Code/Source/svFSI/consts.cpp b/Code/Source/svFSI/consts.cpp
index 59226ed3..d720ab60 100644
--- a/Code/Source/svFSI/consts.cpp
+++ b/Code/Source/svFSI/consts.cpp
@@ -75,6 +75,14 @@ const std::map constitutive_model_name_to_typ
};
+// Map for contact model string name to ContacteModelType.
+const std::map contact_model_name_to_type =
+{
+ {"penalty", ContactModelType::cntctM_penalty},
+ {"potential", ContactModelType::cntctM_potential},
+};
+
+
//------------------------------------
// fluid_viscosity_model_name_to_type
//------------------------------------
diff --git a/Code/Source/svFSI/consts.h b/Code/Source/svFSI/consts.h
index 8f161764..259eb193 100644
--- a/Code/Source/svFSI/consts.h
+++ b/Code/Source/svFSI/consts.h
@@ -190,6 +190,8 @@ enum class ConstitutiveModelType
stIso_lin = 606,
stIso_Gucci = 607,
stIso_HO = 608,
+ stIso_HO_ma = 610,
+ stIso_LS = 611,
stVol_NA = 650,
stVol_Quad = 651,
stVol_ST91 = 652,
@@ -199,6 +201,20 @@ enum class ConstitutiveModelType
// Map for constitutive_model string to ConstitutiveModelType.
extern const std::map constitutive_model_name_to_type;
+//------------------
+// ContactModelType
+//------------------
+//
+enum class ContactModelType
+{
+ cntctM_NA = 800,
+ cntctM_penalty = 801,
+ cntctM_potential = 802
+};
+
+// Map for model type string to ContactModelType.
+extern const std::map contact_model_name_to_type;
+
// Differenty type of coupling for cplBC.
//
// INTEGER(KIND=IKIND), PARAMETER :: cplBC_NA = 400, cplBC_I = 401,
@@ -340,6 +356,9 @@ enum class OutputType
outGrp_strain = 520,
outGrp_divV = 521,
outGrp_Visc = 522,
+ outGrp_fS = 523,
+ outGrp_C = 524,
+ outGrp_I1 = 525,
out_velocity = 599,
out_pressure = 598,
@@ -365,7 +384,10 @@ enum class OutputType
out_defGrad = 578,
out_strain = 577,
out_divergence = 576,
- out_viscosity = 575
+ out_viscosity = 575,
+ out_fibStrn = 574,
+ out_CGstrain = 573,
+ out_CGInv1 = 572
};
//---------------------
diff --git a/Code/Source/svFSI/contact.cpp b/Code/Source/svFSI/contact.cpp
index 171bb641..a36fc34e 100644
--- a/Code/Source/svFSI/contact.cpp
+++ b/Code/Source/svFSI/contact.cpp
@@ -17,9 +17,9 @@ namespace contact {
// contact_forces
//----------------
//
-// [TODO:DaveP] this is not fully implemented.
+// Reproduces Fortran CONSTRUCT_CONTACTPNLTY.
//
-void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
+void construct_contact_pnlty(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
{
using namespace consts;
@@ -31,6 +31,12 @@ void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
const auto& eq = com_mod.eq[cEq];
const auto& cntctM = com_mod.cntctM;
+ #define n_debug_construct_contact_pnlty
+ #ifdef debug_construct_contact_pnlty
+ DebugMsg dmsg(__func__, com_mod.cm.idcm());
+ dmsg.banner();
+ #endif
+
if (eq.phys != EquationType::phys_shell) {
return;
}
@@ -40,6 +46,11 @@ void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
int k = j + 1;
double kl = cntctM.k;
double hl = cntctM.h;
+ #ifdef debug_construct_contact_pnlty
+ //dmsg << "kl: " << kl;
+ //dmsg << "hl: " << hl;
+ //dmsg << "cntctM.c: " << cntctM.c;
+ #endif
// Compute normal vectors at each node in the current configuration
//
@@ -69,7 +80,6 @@ void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
Vector nV1(nsd);
nn::gnns(nsd, eNoN, Nx, xl, nV1, gCov, gCnv);
- //CALL GNNS(eNoN, Nx, xl, nV1, gCov, gCnv)
double Jac = sqrt(utils::norm(nV1));
nV1 = nV1 / Jac;
@@ -83,7 +93,6 @@ void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
for (int i = 0; i < nsd; i++) {
sF(i,Ac) = sF(i,Ac) + w*N(a)*nV1(i);
}
- //sF(:,Ac) = sF(:,Ac) + w*N(a)*nV1(:)
}
}
}
@@ -91,43 +100,35 @@ void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
all_fun::commu(com_mod, sF);
all_fun::commu(com_mod, sA);
- //CALL COMMU(sF)
- //CALL COMMU(sA)
for (int Ac = 0; Ac < tnNo; Ac++) {
if (!utils::is_zero(sA(Ac))) {
for (int i = 0; i < nsd; i++) {
sF(i,Ac) = sF(i,Ac) / sA(Ac);
}
- //sF(:,Ac) = sF(:,Ac)/sA(Ac)
}
double Jac = sqrt(sF.col(Ac) * sF.col(Ac));
- //Jac = sqrt(SUM(sF(:,Ac)**2))
if (!utils::is_zero(Jac)) {
for (int i = 0; i < nsd; i++) {
sF(i,Ac) = sF(i,Ac) / Jac;
}
- //sF(:,Ac) = sF(:,Ac) / Jac
}
}
// Create a bounding box around possible region of contact and bin
// the box with neighboring nodes
//
- // [TODO:DaveP] I'm not sure what to do here, what this code
- // is actually doing. Probably just move it to a subroutine.
- //
int maxNnb = 15;
-// 101 maxNnb = maxNnb + 5
-// IF (ALLOCATED(bBox)) DEALLOCATE(bBox)
-// ALLOCATE(bBox(maxNnb,tnNo))
-// bBox = 0
+ label_101: maxNnb = maxNnb + 5;
+
Array bBox(maxNnb,tnNo);
+ bBox = -1;
for (int iM = 0; iM < com_mod.nMsh; iM++) {
auto& msh = com_mod.msh[iM];
+
if (!msh.lShl) {
continue;
}
@@ -138,7 +139,7 @@ void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
int Ac = msh.gN(a);
x1(0) = com_mod.x(0,Ac) + Dg(i,Ac);
x1(1) = com_mod.x(1,Ac) + Dg(j,Ac);
- x1(1) = com_mod.x(2,Ac) + Dg(k,Ac);
+ x1(2) = com_mod.x(2,Ac) + Dg(k,Ac);
// Box limits for each node
xmin = x1 - cntctM.c;
@@ -162,8 +163,11 @@ void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
if ((x2(0) >= xmin(0)) && (x2(0) <= xmax(0)) && (x2(1) >= xmin(1)) && (x2(1) <= xmax(1)) &&
(x2(2) >= xmin(2)) && (x2(2) <= xmax(2))) {
- for (int l = 0; l < maxNnb; l++) {
- if (bBox(l,Ac) == 0) {
+
+ int l = 0;
+
+ for (int i = 0; i < maxNnb; i++, l++) {
+ if (bBox(l,Ac) == -1) {
bBox(l,Ac) = Bc;
break;
}
@@ -174,15 +178,16 @@ void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
if (Bc == bBox(l,Ac)) {
break;
}
- //if (bBox(maxNnb,Ac) .NE. 0) GOTO 101
- for (int m = maxNnb; m >= l+1; m--) {
+ if (bBox(maxNnb-1,Ac) != -1) goto label_101;
+
+ for (int m = maxNnb-1; m >= l; m--) {
bBox(m,Ac) = bBox(m-1,Ac);
}
bBox(l,Ac) = Bc;
break;
}
- //if (l .GT. maxNnb) GOTO 101
+ if (l > maxNnb-1) goto label_101;
}
} // b
} // jM
@@ -193,28 +198,28 @@ void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
// corresponding penalty forces assembled to the residue
//
Array lR(dof,tnNo);
- Vector incNd(tnNo);
+ Vector incNd(tnNo);
Vector x1(nsd), x2(nsd);
for (int Ac = 0; Ac < tnNo; Ac++) {
- if (bBox(0,Ac) == 0) {
+ if (bBox(0,Ac) == -1) {
continue;
}
x1(0) = com_mod.x(0,Ac) + Dg(i,Ac);
x1(1) = com_mod.x(1,Ac) + Dg(j,Ac);
x1(2) = com_mod.x(2,Ac) + Dg(k,Ac);
- auto nV1 = sF.col(Ac);
+ auto nV1 = sF.rcol(Ac);
int nNb = 0;
for (int a = 0; a < maxNnb; a++) {
int Bc = bBox(a,Ac);
- if (Bc == 0) {
+ if (Bc == -1) {
continue;
}
- x2(0) = com_mod.x(0,Bc) + Dg(i,Bc);
- x2(1) = com_mod.x(1,Bc) + Dg(j,Bc);
- x2(2) = com_mod.x(2,Bc) + Dg(k,Bc);
- auto nV2 = sF.col(Bc);
+ x2(0) = com_mod.x(0,Bc) + Dg(i,Bc);
+ x2(1) = com_mod.x(1,Bc) + Dg(j,Bc);
+ x2(2) = com_mod.x(2,Bc) + Dg(k,Bc);
+ auto nV2 = sF.rcol(Bc);
auto x12 = x1 - x2;
double c = sqrt(utils::norm(x12));
@@ -224,12 +229,12 @@ void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
double d = utils::norm(x12, nV2);
bool flag = false;
double pk{0.0};
-
+
if (d >= -cntctM.h && d < 0.0) {
- pk = 0.5*kl/hl * pow(d + hl,2.0);
+ pk = 0.5 * (kl / hl) * pow(d+hl, 2.0);
flag = true;
} else if (d >= 0.0) {
- pk = 0.5*kl * (hl + d);
+ pk = 0.5 * kl * (hl + d);
flag = true;
} else {
pk = 0.0;
@@ -241,17 +246,26 @@ void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
for (int i = 0; i < nsd; i++) {
lR(i,Ac) = lR(i,Ac) - pk*nV1(i);
}
- //lR(1:nsd,Ac) = lR(1:nsd,Ac) - pk*nV1(:)
+ #ifdef debug_construct_contact_pnlty
+ dmsg << " " << " ";
+ dmsg << "Ac: " << Ac+1;
+ dmsg << "Bc: " << Bc+1;
+ dmsg << "nV1: " << nV1;
+ dmsg << "nV2: " << nV2;
+ dmsg << "pk: " << pk;
+ dmsg << "x12: " << x12;
+ dmsg << "d: " << d;
+ #endif
}
}
}
- if (nNb != 0) {
- for (int i = 0; i < dof; i++) {
- lR(i,Ac) = lR(i,Ac) / static_cast(nNb);
- }
- //lR(:,Ac) = lR(:,Ac) / REAL(nNb, KIND=RKIND)
+ if (nNb != 0) {
+ for (int i = 0; i < dof; i++) {
+ lR(i,Ac) = lR(i,Ac) / static_cast(nNb);
+ }
}
+
}
// Return if no penalty forces are to be added
@@ -260,6 +274,7 @@ void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
}
// Global assembly
+ //
for (int Ac = 0; Ac < tnNo; Ac++) {
if (incNd(Ac) == 0) {
continue;
@@ -267,8 +282,8 @@ void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg)
for (int i = 0; i < dof; i++) {
com_mod.R(i,Ac) = com_mod.R(i,Ac) + lR(i,Ac);
}
- //R(:,Ac) = R(:,Ac) + lR(:,Ac)
}
+
}
};
diff --git a/Code/Source/svFSI/contact.h b/Code/Source/svFSI/contact.h
index df0fa067..db62ab27 100644
--- a/Code/Source/svFSI/contact.h
+++ b/Code/Source/svFSI/contact.h
@@ -5,7 +5,7 @@
namespace contact {
-void contact_forces(ComMod& com_mod, CmMod& cm_mod, const Array& Dg);
+void construct_contact_pnlty(ComMod& com_mod, CmMod& cm_mod, const Array& Dg);
};
diff --git a/Code/Source/svFSI/distribute.cpp b/Code/Source/svFSI/distribute.cpp
index 2166f8c8..0815d19b 100644
--- a/Code/Source/svFSI/distribute.cpp
+++ b/Code/Source/svFSI/distribute.cpp
@@ -312,7 +312,7 @@ void distribute(Simulation* simulation)
cm.bcast(cm_mod, &com_mod.iCntct);
if (com_mod.iCntct) {
- cm.bcast(cm_mod, &com_mod.cntctM.cType);
+ cm.bcast_enum(cm_mod, &com_mod.cntctM.cType);
cm.bcast(cm_mod, &com_mod.cntctM.k);
cm.bcast(cm_mod, &com_mod.cntctM.c);
cm.bcast(cm_mod, &com_mod.cntctM.h);
@@ -1151,6 +1151,11 @@ void dist_mat_consts(const ComMod& com_mod, const CmMod& cm_mod, const cmType& c
cm.bcast(cm_mod, &lStM.afs);
cm.bcast(cm_mod, &lStM.bfs);
cm.bcast(cm_mod, &lStM.kap);
+ cm.bcast(cm_mod, &lStM.khs);
+ cm.bcast(cm_mod, &lStM.a0);
+ cm.bcast(cm_mod, &lStM.b1);
+ cm.bcast(cm_mod, &lStM.b2);
+ cm.bcast(cm_mod, &lStM.mu0);
// Distribute fiber stress
cm.bcast(cm_mod, &lStM.Tf.fType);
diff --git a/Code/Source/svFSI/eq_assem.cpp b/Code/Source/svFSI/eq_assem.cpp
index d0f5fb69..47b6fb7d 100644
--- a/Code/Source/svFSI/eq_assem.cpp
+++ b/Code/Source/svFSI/eq_assem.cpp
@@ -16,6 +16,7 @@
#include "heats.h"
#include "l_elas.h"
#include "mesh.h"
+#include "shells.h"
#include "stokes.h"
#include "sv_struct.h"
#include "ustruct.h"
@@ -352,7 +353,7 @@ void global_eq_assem(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const
break;
case EquationType::phys_shell:
- throw std::runtime_error("[global_eq_assem] phys_shell not implemented.");
+ shells::construct_shell(com_mod, lM, Ag, Yg, Dg);
break;
case EquationType::phys_FSI:
diff --git a/Code/Source/svFSI/fft.cpp b/Code/Source/svFSI/fft.cpp
index d1ecce0a..764e177e 100644
--- a/Code/Source/svFSI/fft.cpp
+++ b/Code/Source/svFSI/fft.cpp
@@ -20,8 +20,9 @@ void fft(const int np, const std::vector>& temporal_values,
#define n_debug_fft
#ifdef debug_fft
- DebugMsg dmsg(__func__, com_mod.cm.idcm());
+ DebugMsg dmsg(__func__, 0);
dmsg.banner();
+ dmsg << "np: " << np;
#endif
Vector t(np);
@@ -62,11 +63,6 @@ void fft(const int np, const std::vector>& temporal_values,
int ns = 0;
Vector ns_array(512);
- double r[1][16];
- for (int n = 0; n < gt.n; n++) {
- r[0][n] = 0.0;
- }
-
for (int n = 0; n < gt.n; n++) {
double tmp = static_cast(n);
gt.r.set_col(n, 0.0);
@@ -80,9 +76,7 @@ void fft(const int np, const std::vector>& temporal_values,
double s = (q(j,i+1) - q(j,i)) / (t[i+1] - t[i]);
if (n == 0) {
gt.r(j,n) = gt.r(j,n) + 0.5*(t[i+1] - t[i]) * (q(j,i+1) + q(j,i));
- r[j][n] = r[j][n] + 0.5*(t[i+1] - t[i]) * (q(j,i+1) + q(j,i));
} else {
- r[j][n] = r[j][n] + s*(cos(kn) - cos(ko));
gt.r(j,n) = gt.r(j,n) + s*(cos(kn) - cos(ko));
gt.i(j,n) = gt.i(j,n) - s*(sin(kn) - sin(ko));
}
diff --git a/Code/Source/svFSI/fluid.cpp b/Code/Source/svFSI/fluid.cpp
index e4ff4a50..2acdcec2 100644
--- a/Code/Source/svFSI/fluid.cpp
+++ b/Code/Source/svFSI/fluid.cpp
@@ -1568,7 +1568,7 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e
ux(1,2) = ux(1,2) + Nwx(1,a)*yl(2,a);
ux(2,2) = ux(2,2) + Nwx(2,a)*yl(2,a);
- uxx(0,0,1) = uxx(0,0,0) + Nwxx(0,a)*yl(0,a);
+ uxx(0,0,0) = uxx(0,0,0) + Nwxx(0,a)*yl(0,a);
uxx(1,0,1) = uxx(1,0,1) + Nwxx(1,a)*yl(0,a);
uxx(2,0,2) = uxx(2,0,2) + Nwxx(2,a)*yl(0,a);
uxx(1,0,0) = uxx(1,0,0) + Nwxx(3,a)*yl(0,a);
diff --git a/Code/Source/svFSI/lhsa.cpp b/Code/Source/svFSI/lhsa.cpp
index bf1b5ad7..899a1a4b 100644
--- a/Code/Source/svFSI/lhsa.cpp
+++ b/Code/Source/svFSI/lhsa.cpp
@@ -177,6 +177,10 @@ void lhsa(Simulation* simulation, int& nnz)
continue;
}
+ if (msh.eType != ElementType::TRI3) {
+ continue;
+ }
+
if (msh.eType == ElementType::NRB) {
continue;
}
@@ -185,7 +189,7 @@ void lhsa(Simulation* simulation, int& nnz)
for (int a = 0; a < 2*msh.eNoN; a++) {
int rowN;
- if (a <= msh.eNoN) {
+ if (a < msh.eNoN) {
rowN = msh.IEN(a,e);
} else {
rowN = msh.eIEN(a-msh.eNoN,e);
@@ -198,7 +202,7 @@ void lhsa(Simulation* simulation, int& nnz)
int colN;
for (int b = 0; b < 2*msh.eNoN; b++) {
- if (b <= msh.eNoN) {
+ if (b < msh.eNoN) {
colN = msh.IEN(b,e);
} else {
colN = msh.eIEN(b-msh.eNoN,e);
diff --git a/Code/Source/svFSI/main.cpp b/Code/Source/svFSI/main.cpp
index accf4faf..90c7c94b 100644
--- a/Code/Source/svFSI/main.cpp
+++ b/Code/Source/svFSI/main.cpp
@@ -328,7 +328,15 @@ void iterate_solution(Simulation* simulation)
// Apply contact model and add its contribution to residue
//
if (com_mod.iCntct) {
- contact::contact_forces(com_mod, cm_mod, Dg);
+ contact::construct_contact_pnlty(com_mod, cm_mod, Dg);
+
+#if 0
+ if (cTS <= 2050) {
+ Array::write_enabled = true;
+ com_mod.R.write("R_"+ std::to_string(cTS));
+ //exit(0);
+ }
+#endif
}
// Synchronize R across processes. Note: that it is important
@@ -648,15 +656,24 @@ int main(int argc, char *argv[])
// Read in the solver commands .xml file.
//
+ #ifdef debug_main
+ dmsg << "Read files " << " ... ";
+ #endif
read_files(simulation, file_name);
// Distribute data to processors.
+ #ifdef debug_main
+ dmsg << "Distribute data to processors " << " ... ";
+ #endif
distribute(simulation);
// Initialize simulation data.
//
Vector init_time(3);
+ #ifdef debug_main
+ dmsg << "Initialize " << " ... ";
+ #endif
initialize(simulation, init_time);
#ifdef debug_main
diff --git a/Code/Source/svFSI/mat_fun.cpp b/Code/Source/svFSI/mat_fun.cpp
index f8a8db88..bfa16349 100644
--- a/Code/Source/svFSI/mat_fun.cpp
+++ b/Code/Source/svFSI/mat_fun.cpp
@@ -649,6 +649,32 @@ double mat_trace(const Array& A, const int nd)
return result;
}
+//-----------------
+// ten_asym_prod12
+//-----------------
+// Create a 4th order tensor from antisymmetric outer product of
+// two matrices
+//
+// Cijkl = Aij*Bkl-Ail*Bjk
+//
+Tensor4
+ten_asym_prod12(const Array& A, const Array& B, const int nd)
+{
+ Tensor4 C(nd,nd,nd,nd);
+
+ int nn = pow(nd,4);
+
+ for (int ii = 0; ii < nn; ii) {
+ int i = t_ind(0,ii);
+ int j = t_ind(1,ii);
+ int k = t_ind(2,ii);
+ int l = t_ind(3,ii);
+ C(i,j,k,l) = 0.5 * ( A(i,j)*B(k,l) - A(i,l)*B(j,k) );
+ }
+
+ return C;
+}
+
//----------
// ten_ddot
//----------
@@ -864,6 +890,39 @@ ten_ids(const int nd)
return A;
}
+//-----------
+// ten_mddot
+//-----------
+// Double dot product of a 4th order tensor and a 2nd order tensor
+//
+// C_ij = (A_ijkl * B_kl)
+//
+Array
+ten_mddot(const Tensor4& A, const Array& B, const int nd)
+{
+ Array C(nd,nd);
+
+ if (nd == 2) {
+ for (int i = 0; i < nd; i++) {
+ for (int j = 0; j < nd; j++) {
+ C(i,j) = A(i,j,0,0)*B(0,0) + A(i,j,0,1)*B(0,1) + A(i,j,1,0)*B(1,0) + A(i,j,1,1)*B(1,1);
+ }
+ }
+
+ } else {
+ for (int i = 0; i < nd; i++) {
+ for (int j = 0; j < nd; j++) {
+ C(i,j) = A(i,j,0,0)*B(0,0) + A(i,j,0,1)*B(0,1) + A(i,j,0,2)*B(0,2) + A(i,j,1,0)*B(1,0) +
+ A(i,j,1,1)*B(1,1) + A(i,j,1,2)*B(1,2) + A(i,j,2,0)*B(2,0) + A(i,j,2,1)*B(2,1) +
+ A(i,j,2,2)*B(2,2);
+ }
+ }
+ }
+
+ return C;
+}
+
+
//---------------
// ten_symm_prod
//---------------
diff --git a/Code/Source/svFSI/mat_fun.h b/Code/Source/svFSI/mat_fun.h
index 6531a9bb..e168b035 100644
--- a/Code/Source/svFSI/mat_fun.h
+++ b/Code/Source/svFSI/mat_fun.h
@@ -34,12 +34,15 @@ namespace mat_fun {
Array mat_symm_prod(const Vector& u, const Vector& v, const int nd);
double mat_trace(const Array& A, const int nd);
+ Tensor4 ten_asym_prod12(const Array& A, const Array& B, const int nd);
Tensor4 ten_ddot(const Tensor4& A, const Tensor4& B, const int nd);
Tensor4 ten_ddot_2412(const Tensor4& A, const Tensor4& B, const int nd);
Tensor4 ten_ddot_3424(const Tensor4& A, const Tensor4& B, const int nd);
Tensor4 ten_dyad_prod(const Array& A, const Array& B, const int nd);
Tensor4 ten_ids(const int nd);
+ Array ten_mddot(const Tensor4& A, const Array& B, const int nd);
+
Tensor4 ten_symm_prod(const Array& A, const Array& B, const int nd);
Tensor4 ten_transpose(const Tensor4& A, const int nd);
diff --git a/Code/Source/svFSI/mat_models.cpp b/Code/Source/svFSI/mat_models.cpp
index c858022f..c43a2927 100644
--- a/Code/Source/svFSI/mat_models.cpp
+++ b/Code/Source/svFSI/mat_models.cpp
@@ -780,6 +780,572 @@ void get_pk2cc_dev(const ComMod& com_mod, const CepMod& cep_mod, const dmnType&
cc_to_voigt(nsd, CC, Dm);
}
+//----------------
+// get_pk2cc_shlc
+//----------------
+// Compute 2nd Piola-Kirchhoff stress and material stiffness tensors
+// for compressible shell elements.
+//
+void get_pk2cc_shlc(const ComMod& com_mod, const dmnType& lDmn, const int nfd, const Array& fNa0,
+ const Array& gg_0, const Array& gg_x, double& g33, Vector& Sml, Array& Dml)
+{
+ // [NOTE] The tolerance here is a bit larger than Fortran.
+ const double ATOL = 1.0e-9;
+ //const double ATOL = 1E-10;
+ const int MAXITR = 20;
+
+ using namespace consts;
+ using namespace mat_fun;
+ using namespace utils;
+
+ #define n_debug_get_pk2cc_shlc
+ #ifdef debug_get_pk2cc_shlc
+ DebugMsg dmsg(__func__, com_mod.cm.idcm());
+ dmsg.banner();
+ #endif
+
+ int nsd = com_mod.nsd;
+ Sml = 0.0;
+ Dml = 0.0;
+
+ // Initialize tensor operations
+ ten_init(3);
+
+ // Some preliminaries
+ auto stM = lDmn.stM;
+ auto kap = stM.Kpen;
+ auto mu = 2.0 * stM.C10;
+ auto f13 = 1.0 / 3.0;
+ auto f23 = 2.0 / 3.0;
+ #ifdef debug_get_pk2cc_shlc
+ dmsg << "kap: " << kap;
+ dmsg << "mu: " << mu;
+ #endif
+
+ // Inverse of metric coefficients in shell continuum
+ auto gi_x = mat_inv(gg_x, 2);
+
+ Array gi_0(3,3);
+ auto gg_0_inv = mat_inv(gg_0, 2);
+ for (int i = 0; i < 2; i++) {
+ for (int j = 0; j < 2; j++) {
+ gi_0(i,j) = gg_0_inv(i,j);
+ }
+ }
+
+ gi_0(2,2) = 1.0;
+
+ // Ratio of inplane Jacobian determinant squared
+ auto Jg2 = mat_det(gg_x, 2) / mat_det(gg_0, 2);
+
+ #ifdef debug_get_pk2cc_shlc
+ dmsg << "gi_0: " << gi_0;
+ dmsg << "gi_x: " << gi_x;
+ dmsg << "Jg2: " << Jg2;
+ #endif
+
+ // Begin Newton iterations to satisfy plane-stress condition.
+ // The objective is to find C33 that satisfies S33 = 0.
+ //
+ int itr = 0;
+ double C33 = 1.0;
+ Array S(3,3);
+ Tensor4 CC(3,3,3,3);
+
+ while (true) {
+ itr = itr + 1;
+ //dmsg << "------- itr: " << itr;
+
+ // Trace (C)
+ auto trC3 = (gg_x(0,0)*gi_0(0,0) + gg_x(0,1)*gi_0(0,1) + gg_x(1,0)*gi_0(1,0) +
+ gg_x(1,1)*gi_0(1,1) + C33)*f13;
+
+ // Jacobian-related quantities
+ auto J2 = Jg2*C33;
+ auto J23 = pow(J2,-f13);
+ //dmsg << "trC3: " << trC3;
+ //dmsg << "J2: " << J2;
+ //dmsg << "J23: " << J23;
+ //dmsg << "f13: " << f13;
+
+ // Inverse of curvilinear Cauchy-Green deformation tensor
+ //
+ Array Ci(3,3);
+ Ci(2,2) = 1.0 / C33;
+
+ for (int i = 0; i < 2; i++) {
+ for (int j = 0; j < 2; j++) {
+ Ci(i,j) = gi_x(i,j);
+ }
+ }
+ //dmsg << "Ci: " << Ci;
+
+ // Contribution from dilational penalty terms to S and CC
+ auto pJ = 0.50 * kap * (J2 - 1.0);
+ auto plJ = kap * J2;
+ #ifdef debug_get_pk2cc_shlc
+ dmsg << "pJ: " << pJ;
+ dmsg << "plJ: " << plJ;
+ dmsg << "J23: " << J23;
+ dmsg << "mu: " << mu;
+ dmsg << "trC3: " << trC3;
+ dmsg << "stM.isoType: " << stM.isoType;
+ #endif
+
+ switch (stM.isoType) {
+
+ case ConstitutiveModelType::stIso_nHook: {
+ // 2nd Piola Kirchhoff stress
+ S = mu*J23*(gi_0 - trC3*Ci) + pJ*Ci;
+
+ // Elasticity tensor
+ CC = (mu*J23*f23*trC3 + plJ)*ten_dyad_prod(Ci, Ci, 3) + (mu*J23*trC3 - pJ)*2.0*ten_symm_prod(Ci, Ci, 3) -
+ f23*mu*J23*(ten_dyad_prod(gi_0, Ci, 3) + ten_dyad_prod(Ci, gi_0, 3));
+ } break;
+
+ case ConstitutiveModelType::stIso_MR: {
+
+ // 2nd Piola Kirchhoff stress
+ auto C1 = stM.C10;
+ auto C2 = stM.C01;
+ auto J43 = pow(J2,-f23);
+
+ auto I2ijkl = ten_dyad_prod(gi_0, gi_0, 3) - ten_symm_prod(gi_0, gi_0, 3);
+ auto I2ij = ten_mddot(I2ijkl, Ci, 3);
+
+ auto Gi4AS = ten_asym_prod12(gi_0, gi_0, 3);
+ auto I2 = mat_ddot(Ci, ten_mddot(Gi4AS, Ci, 3), 3);
+ auto Cikl = -1.0 * ten_symm_prod(Ci, Ci, 3);
+
+ S = C1*J23*(gi_0 - trC3*Ci) + pJ*Ci + C2*J43*(I2ij - f23*I2*Ci);
+
+ // Elasticity tensor
+ CC = (C1*J23*f23*trC3 + plJ)*ten_dyad_prod(Ci, Ci, 3) + (C1*J23*trC3 - pJ)*2.0*ten_symm_prod(Ci, Ci, 3) -
+ f23*C1*J23*(ten_dyad_prod(gi_0, Ci, 3) + ten_dyad_prod(Ci, gi_0, 3));
+
+ CC += 2.0 * f23 * C2 * J43 * ( ten_dyad_prod((f23*I2*Ci-I2ij), Ci, 3) - I2*Cikl -
+ ten_dyad_prod(Ci, I2ij, 3) + I2ijkl);
+ } break;
+
+ case ConstitutiveModelType::stIso_HO_ma: {
+
+ if (nfd != 2) {
+ throw std::runtime_error("Min fiber directions not defined for Holzapfel material model (1)");
+ }
+
+ Array3 fl(2,2,nfd);
+
+ for (int iFn = 0; iFn < nfd; iFn++) {
+ fl(0,0,iFn) = fNa0(0,iFn)*fNa0(0,iFn);
+ fl(0,1,iFn) = fNa0(0,iFn)*fNa0(1,iFn);
+ fl(1,0,iFn) = fNa0(1,iFn)*fNa0(0,iFn);
+ fl(1,1,iFn) = fNa0(1,iFn)*fNa0(1,iFn);
+ }
+
+ // Compute fiber-based invariants
+ //
+ auto Inv4 = gg_x(0,0)*fl(0,0,0) + gg_x(0,1)*fl(0,1,0) + gg_x(1,0)*fl(1,0,0) + gg_x(1,1)*fl(1,1,0);
+ auto Inv6 = gg_x(0,0)*fl(0,0,1) + gg_x(0,1)*fl(0,1,1) + gg_x(1,0)*fl(1,0,1) + gg_x(1,1)*fl(1,1,1);
+ auto Inv8 = gg_x(0,0)*fNa0(0,0)*fNa0(0,1) + gg_x(0,1)*fNa0(0,0)*fNa0(1,1) +
+ gg_x(1,0)*fNa0(1,0)*fNa0(0,1) + gg_x(1,1)*fNa0(1,0)*fNa0(1,1);
+
+ auto Eff = Inv4 - 1.0;
+ auto Ess = Inv6 - 1.0;
+ auto Efs = Inv8;
+
+ // Smoothed heaviside function
+ auto c4f = 1.0 / (1.0 + exp(-stM.khs*Eff));
+ auto c4s = 1.0 / (1.0 + exp(-stM.khs*Ess));
+
+ // Approx. derivative of smoothed heaviside function
+ auto dc4f = 0.250 * stM.khs * exp(-stM.khs*fabs(Eff));
+ auto dc4s = 0.250 * stM.khs * exp(-stM.khs*fabs(Ess));
+
+ // Add isochoric stress and stiffness contribution
+ //
+ // EI1 = I1 + Jg2i - 3.0
+ //
+ auto d1 = stM.a*J23*exp(2.0*stM.b*(trC3*J23 - 1.0));
+ auto SN = (gi_0 - trC3*Ci);
+
+ S = d1*SN + pJ*Ci;
+
+ CC = f23*trC3*ten_dyad_prod(Ci, Ci, 3)
+ + trC3*2.0*ten_symm_prod(Ci, Ci, 3)
+ - f23*(ten_dyad_prod(gi_0, Ci, 3)
+ + ten_dyad_prod(Ci, gi_0, 3))
+ + 2.0*stM.b*J23*ten_dyad_prod(SN, SN, 3);
+
+ CC = d1*CC + plJ*ten_dyad_prod(Ci, Ci, 3) - pJ*2.0*ten_symm_prod(Ci, Ci, 3);
+
+ // Anisotropic part
+ // Fiber sheet
+ //
+ Array Hfs(3,3);
+ auto hfs_sym_prod = mat_symm_prod(fNa0.col(0), fNa0.col(1), 2);
+
+ for (int i = 0; i < 2; i++) {
+ for (int j = 0; j < 2; j++) {
+ Hfs(i,j) = hfs_sym_prod(i,j);;
+ }
+ }
+
+ auto g1 = 2.0*stM.afs*exp(stM.bfs*Efs*Efs);
+ S += g1*Efs*Hfs;
+
+ auto g2 = g1*2.0*(1.0 + 2.0*Efs*Efs);
+ CC += g2*ten_dyad_prod(Hfs, Hfs, 3);
+
+ // Fiber
+ Array flM(3,3);
+
+ if (Eff > 0.0) {
+ for (int i = 0; i < 2; i++) {
+ for (int j = 0; j < 2; j++) {
+ flM(i,j) = fl(i,j,0);
+ }
+ }
+
+ // S = S + 2.0*stM.aff*Eff*flM
+ // CC = CC + 4.0*stM.aff*ten_dyad_prod(flM, flM, 2)
+ g1 = 2.0*stM.aff*exp(stM.bff*Eff*Eff);
+ S += g1 * Eff * flM;
+ g2 = g1 * 2.0 * (1.0 + 2.0 * stM.bff * Eff * Eff);
+ CC += g2 * ten_dyad_prod(flM, flM, 3);
+ }
+
+ // Sheet
+ //
+ if (Ess > 0.0) {
+ Array flM(3,3);
+ auto flM_1 = fl.rslice(1);
+ for (int i = 0; i < 2; i++) {
+ for (int j = 0; j < 2; j++) {
+ flM(i,j) = flM_1(i,j);
+ }
+ }
+
+ auto g1 = 2.0 * stM.ass * exp(stM.bss*Ess*Ess);
+ S += g1*Ess*flM;
+ auto g2 = g1 * 2.0 * (1.0 + 2.0 * stM.bss * Ess * Ess);
+ CC += g2*ten_dyad_prod(flM, flM, 3);
+ }
+ } break;
+
+ case ConstitutiveModelType::stIso_LS: {
+ Array3 fl(2,2,nfd);
+
+ for (int iFn = 0; iFn < nfd; iFn++) {
+ fl(0,0,iFn) = fNa0(0,iFn)*fNa0(0,iFn);
+ fl(0,1,iFn) = fNa0(0,iFn)*fNa0(1,iFn);
+ fl(1,0,iFn) = fNa0(1,iFn)*fNa0(0,iFn);
+ fl(1,1,iFn) = fNa0(1,iFn)*fNa0(1,iFn);
+ }
+
+ // Compute fiber-based invariants
+
+ auto Inv4 = gg_x(0,0)*fl(0,0,0) + gg_x(0,1)*fl(0,1,0) + gg_x(1,0)*fl(1,0,0) + gg_x(1,1)*fl(1,1,0);
+ auto Eff = Inv4 - 1.0;
+
+ // Isotropic contribution
+ //
+ auto d1 = 2.0*stM.a0*stM.mu0*stM.b1*J23 * exp(2.0*stM.b1*(trC3*J23 - 1.0));
+ auto SN = (gi_0 - trC3*Ci);
+ auto CCN = f23*trC3*ten_dyad_prod(Ci, Ci, 3) + trC3*2.0*ten_symm_prod(Ci, Ci, 3) -
+ f23*(ten_dyad_prod(gi_0, Ci, 3) + ten_dyad_prod(Ci, gi_0, 3));
+
+ S = (stM.a + d1*(2.0*(trC3*J23 - 1.0)))*SN + pJ*Ci;
+
+ CC = (1.0 + 18.0*stM.b1*(trC3*J23 - 1.0) * (trC3*J23 - 1.0))*J23*ten_dyad_prod(SN, SN, 3) +
+ 3.0*(trC3*J23 - 1.0)*CCN;
+ CC = stM.a*CCN + 2.0*d1*CC + plJ*ten_dyad_prod(Ci, Ci, 3) - pJ*2.0*ten_symm_prod(Ci, Ci, 3);
+
+ // Anisotropic contribution
+ Array flM(3,3);
+
+ if (Eff > 0.0) {
+ auto flM_0 = fl.rslice(0);
+ for (int i = 0; i < 2; i++) {
+ for (int j = 0; j < 2; j++) {
+ flM(i,j) = flM_0(i,j);
+ }
+ }
+
+ S += 2.0*stM.a0*(1.0 - stM.mu0) *stM.b2*Eff * exp(stM.b2*Eff*Eff)*flM;
+
+ CC += 4.0*stM.a0*(1.0 - stM.mu0) * stM.b2 * (1.0 + 2.0*stM.b2*Eff*Eff) *
+ exp(stM.b2*Eff*Eff) * ten_dyad_prod(flM, flM, 3);
+ }
+
+ } break;
+
+ default:
+ //err = "Undefined material constitutive model"
+ break;
+
+ }
+
+ if (fabs(S(2,2)) <= ATOL) {
+ break;
+ }
+
+ if (itr > MAXITR) {
+ std::cerr << "[get_pk2cc_shlc] Failed to converge plane-stress condition." << std::endl;
+ //exit(0);
+ break;
+ }
+
+ C33 = C33 - (2.0 * S(2,2) / CC(2,2,2,2));
+ //dmsg << "1: C33: " << C33;
+ //dmsg << "CC(3,3,3,3): " << CC(2,2,2,2);
+ //exit(0);
+ }
+
+ g33 = C33;
+
+ // Statically condense CC
+ //
+ for (int i = 0; i < 2; i++) {
+ for (int j = 0; j < 2; j++) {
+ for (int k = 0; k < 2; k++) {
+ for (int l = 0; l < 2; l++) {
+ C33 = CC(i,j,2,2)*CC(2,2,k,l) / CC(2,2,2,2);
+ CC(i,j,k,l) = CC(i,j,k,l) - C33;
+ }
+ }
+ }
+ }
+
+ g33 = C33;
+
+ //dmsg << "2: C33: " << C33;
+ //exit(0);
+
+ // Convert the in-plane components to Voigt notation
+ Sml(0) = S(0,0);
+ Sml(1) = S(1,1);
+ Sml(2) = S(0,1);
+
+ Dml(0,0) = CC(0,0,0,0);
+ Dml(0,1) = CC(0,0,1,1);
+ Dml(0,2) = CC(0,0,0,1);
+
+ Dml(1,1) = CC(1,1,1,1);
+ Dml(1,2) = CC(1,1,0,1);
+
+ Dml(2,2) = CC(0,1,0,1);
+
+ Dml(1,0) = Dml(0,1);
+ Dml(2,0) = Dml(0,2);
+ Dml(2,1) = Dml(1,2);
+}
+
+//----------------
+// get_pk2cc_shli
+//----------------
+// Compute 2nd Piola-Kirchhoff stress and material stiffness tensors
+// for incompressible shell elements
+//
+// Reproduces Fortran GETPK2CC_SHLi
+//
+void get_pk2cc_shli(const ComMod& com_mod, const dmnType& lDmn, const int nfd, const Array& fNa0,
+ const Array& gg_0, const Array& gg_x, double& g33, Vector& Sml, Array& Dml)
+{
+ using namespace consts;
+ using namespace mat_fun;
+ using namespace utils;
+
+ #define n_debug_get_pk2cc_shli
+ #ifdef debug_get_pk2cc_shli
+ DebugMsg dmsg(__func__, com_mod.cm.idcm());
+ dmsg.banner();
+ #endif
+
+ int nsd = com_mod.nsd;
+ Sml = 0.0;
+ Dml = 0.0;
+
+ // Some preliminaries
+ auto stM = lDmn.stM;
+
+ // Inverse of metric coefficients in shell continuum
+ auto gi_0 = mat_inv(gg_0,2);
+ auto gi_x = mat_inv(gg_x,2);
+
+ // Ratio of inplane Jacobian determinants
+ auto Jg2i = mat_det(gg_x,2);
+
+ if (is_zero(Jg2i)) {
+ throw std::runtime_error(" Divide by zero in-plane Jacobian determinant.");
+ }
+
+ Jg2i = mat_det(gg_0,2) / Jg2i;
+
+ double I1 = 0.0;
+
+ for (int a = 0; a < 2; a++) {
+ for (int b = 0; b < 2; b++) {
+ I1 = I1 + gi_0(a,b) * gg_x(a,b);
+ }
+ }
+
+ Tensor4 CC(nsd,nsd,nsd,nsd);
+ Array S;
+ Array3 fl(2,2,nfd);
+
+ switch (stM.isoType) {
+
+ case ConstitutiveModelType::stIso_nHook: {
+ auto mu = 2.0 * stM.C10;
+ S = mu*(gi_0 - Jg2i*gi_x);
+ auto CC = 2.0*mu*Jg2i*(ten_dyad_prod(gi_x, gi_x,1) + ten_symm_prod(gi_x, gi_x,1));
+ } break;
+
+ case ConstitutiveModelType::stIso_MR: {
+ auto SN = (gi_0 - Jg2i*gi_x);
+ auto CCN = 2.0*Jg2i*(ten_dyad_prod(gi_x, gi_x,1) + ten_symm_prod(gi_x, gi_x,1));
+ S = stM.C10*SN + stM.C01*Jg2i* (gi_0 - I1*gi_x) + stM.C01/Jg2i*gi_x;
+
+ CC = (stM.C10 + stM.C01*I1) * CCN - 2.0*stM.C01 * Jg2i * (ten_dyad_prod(gi_0, gi_x,1) +
+ ten_dyad_prod(gi_x, gi_0,1)) + 2.0*stM.C01 / Jg2i *(ten_dyad_prod(gi_x, gi_0,1) -
+ ten_symm_prod(gi_x, gi_x,1));
+ } break;
+
+ // HO (Holzapfel-Ogden) model for myocardium with full invariants
+ // for the anisotropy terms (modified-anisotropy)
+ //
+ case ConstitutiveModelType::stIso_HO_ma: {
+
+ if (nfd != 2) {
+ throw std::runtime_error("Min fiber directions not defined for Holzapfel material model (1)");
+ }
+
+ for (int iFn = 0; iFn < nfd; iFn++) {
+ fl(0,0,iFn) = fNa0(0,iFn)*fNa0(0,iFn);
+ fl(0,1,iFn) = fNa0(0,iFn)*fNa0(1,iFn);
+ fl(1,0,iFn) = fNa0(1,iFn)*fNa0(0,iFn);
+ fl(1,1,iFn) = fNa0(1,iFn)*fNa0(1,iFn);
+ }
+
+ // Compute fiber-based invariants
+ auto Inv4 = gg_x(0,0)*fl(0,0,0) + gg_x(0,1)*fl(0,1,0) + gg_x(1,0)*fl(1,0,0) + gg_x(1,1)*fl(1,1,0);
+ auto Inv6 = gg_x(0,0)*fl(0,0,1) + gg_x(0,1)*fl(0,1,1) + gg_x(1,0)*fl(1,0,1) + gg_x(1,1)*fl(1,1,1);
+ auto Inv8 = gg_x(0,0)*fNa0(0,0)*fNa0(0,1) + gg_x(0,1)*fNa0(0,0)*fNa0(1,1) +
+ gg_x(1,0)*fNa0(1,0)*fNa0(0,1) + gg_x(1,1)*fNa0(1,0)*fNa0(1,1);
+
+ auto Eff = Inv4 - 1.0;
+ auto Ess = Inv6 - 1.0;
+ auto Efs = Inv8;
+
+ // Smoothed heaviside function
+ auto c4f = 1.0 / (1.0 + exp(-stM.khs*Eff));
+ auto c4s = 1.0 / (1.0 + exp(-stM.khs*Ess));
+
+ // Approx. derivative of smoothed heaviside function
+ auto dc4f = 0.250*stM.khs*exp(-stM.khs*fabs(Eff));
+ auto dc4s = 0.250*stM.khs*exp(-stM.khs*fabs(Ess));
+
+ // Add isochoric stress and stiffness contribution
+ //
+ auto EI1 = I1 + Jg2i - 3.0;
+ auto SN = (gi_0 - Jg2i*gi_x);
+ auto CCN = 2.0*Jg2i*(ten_dyad_prod(gi_x, gi_x,1) + ten_symm_prod(gi_x, gi_x,1));
+
+ auto d1 = stM.a*exp(stM.b*EI1);
+ S = d1*SN;
+ CC = d1*(CCN + 2.0*stM.b*ten_dyad_prod(SN, SN,1));
+
+ // Anisotropic part
+ // Fiber sheet
+ auto Hfs = mat_symm_prod(fNa0.col(0), fNa0.col(0),1);
+ auto g1 = 2.0*stM.afs*exp(stM.bfs*Efs*Efs);
+ S += g1*Efs*Hfs;
+
+ auto g2 = g1*2.0*(1.0 + 2.0*Efs*Efs);
+ CC += CC + g2*ten_dyad_prod(Hfs, Hfs,1);
+
+ // Fiber
+ if (Eff > 0.0) {
+ auto flM = fl.slice(0);
+ // S = S + 2.0*stM.aff*Eff*flM
+ // CC = CC + 4.0*stM.aff*ten_dyad_prod(flM, flM,1)
+ g1 = 2.0*stM.aff*exp(stM.bff*Eff*Eff);
+ S += g1*Eff*flM;
+ g2 = g1*2.0*(1.0 + 2.0*stM.bff*Eff*Eff);
+ CC += g2*ten_dyad_prod(flM, flM,1);
+ }
+
+ // Sheet
+ if (Ess > 0.0) {
+ auto flM = fl.slice(0);
+ g1 = 2.0*stM.ass*exp(stM.bss*Ess*Ess);
+ S += g1*Ess*flM;
+ g2 = g1*2.0*(1.0 + 2.0*stM.bss*Ess*Ess);
+ CC += g2*ten_dyad_prod(flM, flM,1);
+ }
+ } break;
+
+ // Lee Sacks model for aorta with full invariants
+ // for the anisotropy terms (modified-anisotropy)
+ //
+ case ConstitutiveModelType::stIso_LS: {
+ auto SN = (gi_0 - Jg2i*gi_x);
+ auto CCN = 2.0*Jg2i*(ten_dyad_prod(gi_x, gi_x,1) + ten_symm_prod(gi_x, gi_x,1));
+
+ for (int iFn = 0; iFn < nfd; iFn++) {
+ fl(0,0,iFn) = fNa0(0,iFn)*fNa0(0,iFn);
+ fl(0,1,iFn) = fNa0(0,iFn)*fNa0(1,iFn);
+ fl(1,0,iFn) = fNa0(1,iFn)*fNa0(0,iFn);
+ fl(1,1,iFn) = fNa0(1,iFn)*fNa0(1,iFn);
+ }
+
+ // Compute fiber-based invariants
+ auto Inv4 = gg_x(0,0)*fl(0,0,0) + gg_x(0,1)*fl(0,1,0) + gg_x(1,0)*fl(1,0,0) + gg_x(1,1)*fl(1,1,0);
+ auto Eff = Inv4 - 1.0;
+
+ // Isotropic contribution
+ auto EI1 = (I1 + Jg2i - 3.0);
+ S = stM.a*SN + 2.0*stM.a0*stM.mu0*stM.b1*EI1 * exp(stM.b1*EI1*EI1) * SN;
+ CC = stM.a*CCN + 4.0*stM.a0*stM.mu0*stM.b1 * exp(stM.b1*EI1*EI1) * ((1.0 + 2.0*stM.b1*EI1*EI1) *
+ ten_dyad_prod(SN, SN,1) + 0.50*EI1*CCN);
+
+ // Anisotropic contribution
+ if (Eff > 0.0) {
+ auto flM = fl.rslice(0);
+ S += 2.0*stM.a0*(1.0 - stM.mu0)*stM.b2*Eff * exp(stM.b2*Eff*Eff)*flM;
+ CC += 4.0*stM.a0*(1.0 - stM.mu0)*stM.b2 * (1.0 + 2.0*stM.b2*Eff*Eff) *
+ exp(stM.b2*Eff*Eff) * ten_dyad_prod(flM, flM,1);
+ }
+ } break;
+
+ default:
+ //err = "Undefined material constitutive model"
+ break;
+ }
+
+ g33 = Jg2i;
+
+ // Convert to Voigt notation
+ Sml(0) = S(0,0);
+ Sml(1) = S(1,1);
+ Sml(2) = S(0,1);
+
+ Dml(0,0) = CC(0,0,0,0);
+ Dml(0,1) = CC(0,0,1,1);
+ Dml(0,2) = CC(0,0,0,1);
+
+ Dml(1,1) = CC(1,1,1,1);
+ Dml(1,2) = CC(1,1,0,1);
+
+ Dml(2,2) = CC(0,1,0,1);
+
+ Dml(1,0) = Dml(0,1);
+ Dml(2,0) = Dml(0,2);
+ Dml(2,1) = Dml(1,2);
+
+}
+
+
//------------
// get_svol_p
//------------
diff --git a/Code/Source/svFSI/mat_models.h b/Code/Source/svFSI/mat_models.h
index 0b09164b..5250f24c 100644
--- a/Code/Source/svFSI/mat_models.h
+++ b/Code/Source/svFSI/mat_models.h
@@ -23,6 +23,12 @@ void get_pk2cc(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn
void get_pk2cc_dev(const ComMod& com_mod, const CepMod& cep_mod, const dmnType& lDmn, const Array& F, const int nfd,
const Array& fl, const double ya, Array& S, Array& Dm, double& Ja);
+void get_pk2cc_shlc(const ComMod& com_mod, const dmnType& lDmn, const int nfd, const Array& fNa0,
+ const Array& gg_0, const Array& gg_x, double& g33, Vector& Sml, Array& Dml);
+
+void get_pk2cc_shli(const ComMod& com_mod, const dmnType& lDmn, const int nfd, const Array& fNa0,
+ const Array& gg_0, const Array& gg_x, double& g33, Vector& Sml, Array& Dml);
+
void get_tau(const ComMod& com_mod, const dmnType& lDmn, const double detF, const double Je, double& tauM, double& tauC);
void get_svol_p(const ComMod& com_mod, const CepMod& cep_mod, const stModelType& stM, const double J,
diff --git a/Code/Source/svFSI/nn.cpp b/Code/Source/svFSI/nn.cpp
index 0959da06..8b1fd56b 100644
--- a/Code/Source/svFSI/nn.cpp
+++ b/Code/Source/svFSI/nn.cpp
@@ -565,8 +565,8 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g,
dmsg << "iM: " << iM+1;
dmsg << "Ec: " << Ec+1;
dmsg << "eNoN: " << eNoN;
- dmsg << "msh.IEN.nrows: " << msh.IEN.nrows_;
- dmsg << "msh.IEN.ncols: " << msh.IEN.ncols_;
+ dmsg << "msh.IEN.nrows: " << msh.IEN.nrows();
+ dmsg << "msh.IEN.ncols: " << msh.IEN.ncols();
#endif
Array lX(nsd,eNoN);
@@ -637,7 +637,7 @@ void gnnb(const ComMod& com_mod, const faceType& lFa, const int e, const int g,
// Compute adjoining mesh element normal
//
- Array xXi(nsd,insd);
+ Array xXi(nsd,nsd-1);
for (int a = 0; a < eNoN; a++) {
for (int i = 0; i < insd; i++) {
diff --git a/Code/Source/svFSI/post.cpp b/Code/Source/svFSI/post.cpp
index acf1580d..f2a4868d 100644
--- a/Code/Source/svFSI/post.cpp
+++ b/Code/Source/svFSI/post.cpp
@@ -8,6 +8,7 @@
#include "mat_fun.h"
#include "mat_models.h"
#include "nn.h"
+#include "shells.h"
#include "utils.h"
#include "vtk_xml.h"
#include
@@ -1197,6 +1198,441 @@ void ppbin2vtk(Simulation* simulation)
MPI_Finalize();
}
+//----------
+// shl_post
+//----------
+// Routine for post processing shell-based quantities
+//
+// Reproduces Fortran SHLPOST.
+//
+void shl_post(Simulation* simulation, const mshType& lM, const int m, Array& res,
+ Vector& resE, const Array& lD, const int iEq, consts::OutputType outGrp)
+{
+ using namespace consts;
+ using namespace mat_fun;
+ using namespace utils;
+
+ #define n_debug_shl_post
+ #ifdef debug_shl_post
+ DebugMsg dmsg(__func__, 0);
+ dmsg.banner();
+ #endif
+
+ auto& com_mod = simulation->com_mod;
+ auto& cm = com_mod.cm;
+ auto& cm_mod = simulation->cm_mod;
+ auto& eq = com_mod.eq[iEq];
+
+ const int nsd = com_mod.nsd;
+ const int tnNo = com_mod.tnNo;
+ const int tDof = com_mod.tDof;
+
+ // [NOTE] Setting gobal variable 'dof'.
+ com_mod.dof = eq.dof;
+
+ int i = eq.s;
+ int j = i + 1;
+ int k = j + 1;
+
+ int nFn = lM.nFn;
+ if (nFn == 0) {
+ nFn = 1;
+ }
+
+ // Set shell dimension := 2
+ int insd = nsd-1;
+
+ // Initialize tensor operations
+ ten_init(insd);
+
+ // Set eNoN (number of nodes per element)
+ //
+ int eNoN = lM.eNoN;
+
+ if (lM.eType == ElementType::TRI3) {
+ eNoN = 2*eNoN;
+ }
+
+ Vector sA(tnNo), sE(lM.nEl), resl(m), N(lM.eNoN);
+ Vector ptr(eNoN);
+ Array sF(m,tnNo), dl(tDof,eNoN), x0(3,eNoN), xc(3,eNoN),
+ fN(3,nFn), fNa0(2,eNoN), Nx(2,lM.eNoN);
+ Array3 Bb(3,3,6);
+
+ // Initialize arrays
+ sA = 0.0;
+ sF = 0.0;
+ sE = 0.0;
+ Bb = 0.0;
+
+ // Compute quantities at the element level and project them to nodes
+ //
+ for (int e = 0; e < lM.nEl; e++) {
+ int cDmn = all_fun::domain(com_mod, lM, iEq, e);
+ auto cPhys = eq.dmn[cDmn].phys;
+ //dmsg << "========== e: " << e+1;
+
+ if (cPhys != EquationType::phys_shell) {
+ continue;
+ }
+ //if (lM.eType .EQ. eType_NRB) CALL NRBNNX(lM, e)
+
+ // Get shell properties
+ double nu = eq.dmn[cDmn].prop.at(PhysicalProperyType::poisson_ratio);
+ double ht = eq.dmn[cDmn].prop.at(PhysicalProperyType::shell_thickness);
+
+ // Check for incompressibility
+ //
+ bool incompFlag = false;
+ if (is_zero(nu-0.50)) {
+ incompFlag = true;
+ }
+ //dmsg << "incompFlag: " << incompFlag;
+
+ // Get the reference configuration and displacement field
+ x0 = 0.0;
+ dl = 0.0;
+
+ for (int a = 0; a < eNoN; a++) {
+ int Ac;
+ if (a < lM.eNoN) {
+ Ac = lM.IEN(a,e);
+ ptr(a) = Ac;
+ } else {
+ int b = a - lM.eNoN;
+ Ac = lM.eIEN(b,e);
+ ptr(a) = Ac;
+ if (Ac == -1) {
+ continue;
+ }
+ }
+
+ for (int i = 0; i < 3; i++) {
+ x0(i,a) = com_mod.x(i,Ac);
+ }
+
+ for (int i = 0; i < tDof; i++) {
+ dl(i,a) = lD(i,Ac);
+ }
+ }
+
+ // Get the current configuration
+ xc = 0.0;
+
+ for (int a = 0; a < eNoN; a++) {
+ xc(0,a) = x0(0,a) + dl(i,a);
+ xc(1,a) = x0(1,a) + dl(j,a);
+ xc(2,a) = x0(2,a) + dl(k,a);
+ }
+
+ // Get fiber directions
+ //
+ fN = 0.0;
+
+ if (lM.fN.size() != 0) {
+ for (int iFn = 0; iFn < nFn; iFn++) {
+ for (int i = 0; i < 3; i++) {
+ fN(i,iFn) = lM.fN(i+nsd*iFn,e);
+ }
+ }
+ }
+
+ // Set number of integration points.
+ // Note: Gauss integration performed for NURBS elements
+ // Not required for constant-strain triangle elements
+ //
+ int nwg;
+
+ if (lM.eType == ElementType::TRI3) {
+ nwg = 1;
+ } else {
+ nwg = lM.nG;
+ }
+
+ // Update shapefunctions for NURBS elements
+ //
+ //if (lM.eType .EQ. eType_NRB) CALL NRBNNX(lM, e)
+
+ double Je = 0.0;
+ resl = 0.0;
+ Vector N;
+ Array Nxx, Nx;
+
+ for (int g = 0; g < nwg; g++) {
+ double w = 0.0;
+ double lam3 = 0.0;
+ double aa_0[2][2]{}, bb_0[2][2]{};
+ double aa_x[2][2]{}, bb_x[2][2]{};
+
+ Array aCov0(3,2), aCnv0(3,2), aCov(3,2), aCnv(3,2);
+ Vector nV0(3), nV(3);
+ //dmsg << "---------- g: " << g;
+
+ // [TODO] This is not fully implemented
+ //
+ if (lM.eType != ElementType::TRI3) {
+ // Set element shape functions and their derivatives
+ if (lM.eType == ElementType::NRB) {
+ //N = lM.N(:,g)
+ //Nx = lM.Nx(:,:,g)
+ //Nxx = lM.Nxx(:,:,g)
+ } else {
+ N = lM.fs[0].N.rcol(g);
+ Nx = lM.fs[0].Nx.rslice(g);
+ Nxx = lM.fs[0].Nxx.rslice(g);
+ }
+
+ // Covariant and contravariant bases (ref. config.)
+ //
+ nn::gnns(nsd, eNoN, Nx, x0, nV0, aCov0, aCnv0);
+ auto Jac0 = sqrt(norm(nV0));
+ nV0 = nV0 / Jac0;
+
+ // Covariant and contravariant bases (spatial config.)
+ nn::gnns(nsd, eNoN, Nx, xc, nV, aCov, aCnv);
+ auto Jac = sqrt(norm(nV));
+ nV = nV/Jac;
+
+ // Second derivatives for curvature coeffs. (ref. config)
+#if 0
+ r0_xx(:,:,:) = 0.0
+ r_xx(:,:,:) = 0.0
+
+ for (int a = 0; a < eNoN; a++) {
+ r0_xx(0,0,:) = r0_xx(0,0,:) + Nxx(0,a)*x0(:,a)
+ r0_xx(1,1,:) = r0_xx(1,1,:) + Nxx(1,a)*x0(:,a)
+ r0_xx(0,1,:) = r0_xx(0,1,:) + Nxx(2,a)*x0(:,a)
+
+ r_xx(0,0,:) = r_xx(0,0,:) + Nxx(0,a)*xc(:,a)
+ r_xx(1,1,:) = r_xx(1,1,:) + Nxx(1,a)*xc(:,a)
+ r_xx(0,1,:) = r_xx(0,1,:) + Nxx(2,a)*xc(:,a)
+ }
+
+ r0_xx(1,0,:) = r0_xx(0,1,:)
+ r_xx(1,0,:) = r_xx(0,1,:)
+
+ // Compute metric tensor (aa) and curvature coefficients(bb)
+ //
+ double aa_0[2][2]{}, bb_0[2][2]{};
+ double aa_x[2][2]{}, bb_x[2][2]{};
+
+ for (int l = 0; l < nsd; l++) {
+ aa_0(0,0) = aa_0(0,0) + aCov0(l,0)*aCov0(l,0)
+ aa_0(0,1) = aa_0(0,1) + aCov0(l,0)*aCov0(l,1)
+ aa_0(1,0) = aa_0(1,0) + aCov0(l,1)*aCov0(l,0)
+ aa_0(1,1) = aa_0(1,1) + aCov0(l,1)*aCov0(l,1)
+
+ aa_x(0,0) = aa_x(0,0) + aCov(l,0)*aCov(l,0)
+ aa_x(0,1) = aa_x(0,1) + aCov(l,0)*aCov(l,1)
+ aa_x(1,0) = aa_x(1,0) + aCov(l,1)*aCov(l,0)
+ aa_x(1,1) = aa_x(1,1) + aCov(l,1)*aCov(l,1)
+
+ bb_0(0,0) = bb_0(0,0) + r0_xx(0,0,l)*nV0(l)
+ bb_0(0,1) = bb_0(0,1) + r0_xx(0,1,l)*nV0(l)
+ bb_0(1,0) = bb_0(1,0) + r0_xx(1,0,l)*nV0(l)
+ bb_0(1,1) = bb_0(1,1) + r0_xx(1,1,l)*nV0(l)
+
+ bb_x(0,0) = bb_x(0,0) + r_xx(0,0,l)*nV(l)
+ bb_x(0,1) = bb_x(0,1) + r_xx(0,1,l)*nV(l)
+ bb_x(1,0) = bb_x(1,0) + r_xx(1,0,l)*nV(l)
+ bb_x(1,1) = bb_x(1,1) + r_xx(1,1,l)*nV(l)
+ }
+
+ // Set weight of the Gauss point
+ w = lM.w(g)*Jac0
+
+#endif
+
+ // for constant strain triangles
+ } else {
+
+ // Set element shape functions and their derivatives
+ N = lM.N.rcol(g);
+ Nx = lM.Nx.rslice(0);
+
+ // Covariant and contravariant bases (ref. config.)
+ //
+ Array tmpX(nsd,lM.eNoN);
+
+ for (int i = 0; i < nsd; i++) {
+ for (int j = 0; j < lM.eNoN; j++) {
+ tmpX(i,j) = x0(i,j);
+ }
+ }
+
+ nn::gnns(nsd, lM.eNoN, Nx, tmpX, nV0, aCov0, aCnv0);
+ auto Jac0 = sqrt(norm(nV0));
+ nV0 = nV0 / Jac0;
+
+ // Covariant and contravariant bases (spatial config.)
+ //
+ for (int i = 0; i < nsd; i++) {
+ for (int j = 0; j < lM.eNoN; j++) {
+ tmpX(i,j) = xc(i,j);
+ }
+ }
+
+ nn::gnns(nsd, lM.eNoN, Nx, tmpX, nV, aCov, aCnv);
+ auto Jac = sqrt(norm(nV));
+ nV = nV / Jac;
+
+ // Compute metric tensor (aa)
+ //
+ for (int l = 0; l < nsd; l++) {
+ aa_0[0][0] = aa_0[0][0] + aCov0(l,0)*aCov0(l,0);
+ aa_0[0][1] = aa_0[0][1] + aCov0(l,0)*aCov0(l,1);
+ aa_0[1][0] = aa_0[1][0] + aCov0(l,1)*aCov0(l,0);
+ aa_0[1][1] = aa_0[1][1] + aCov0(l,1)*aCov0(l,1);
+
+ aa_x[0][0] = aa_x[0][0] + aCov(l,0)*aCov(l,0);
+ aa_x[0][1] = aa_x[0][1] + aCov(l,0)*aCov(l,1);
+ aa_x[1][0] = aa_x[1][0] + aCov(l,1)*aCov(l,0);
+ aa_x[1][1] = aa_x[1][1] + aCov(l,1)*aCov(l,1);
+ }
+
+ shells::shell_bend_cst(com_mod, lM, e, ptr, x0, xc, bb_0, bb_x, Bb, false);
+
+ // Set weight of the Gauss point
+ w = Jac0*0.50;
+ }
+
+ // Compute fiber direction in curvature coordinates
+ //
+ fNa0 = 0.0;
+
+ for (int iFn = 0; iFn < nFn; iFn++) {
+ for (int l = 0; l < nsd; l++) {
+ fNa0(0,iFn) = fNa0(0,iFn) + fN(l,iFn)*aCnv0(l,0);
+ fNa0(1,iFn) = fNa0(1,iFn) + fN(l,iFn)*aCnv0(l,1);
+ }
+ }
+
+ // Compute stress resultants and lambda3 (integrated through
+ // the shell thickness)
+ //
+ Array Sm(3,2);
+ Array3 Dm(3,3,3);
+ shells::shl_strs_res(com_mod, eq.dmn[cDmn], nFn, fNa0, aa_0, aa_x, bb_0, bb_x, lam3, Sm, Dm);
+
+ // Shell in-plane deformation gradient tensor
+ //
+ auto F = mat_dyad_prod(aCov.rcol(0), aCnv0.rcol(0), 3) +
+ mat_dyad_prod(aCov.rcol(1), aCnv0.rcol(1), 3);
+
+ // D deformation gradient tensor in shell continuum
+ auto F3d = F + lam3*mat_dyad_prod(nV, nV0, 3);
+ auto detF = mat_det(F3d, nsd);
+ Je = Je + w;
+ auto Im = mat_id(nsd);
+
+ switch (outGrp) {
+ //dmsg << "outGrp: " << outGrp;
+ case OutputType::outGrp_J: {
+ // Jacobian := determinant of deformation gradient tensor
+ resl(0) = detF;
+ sE(e) = sE(e) + w*detF;
+ } break;
+
+ case OutputType::outGrp_F:
+ // 3D deformation gradient tensor (F)
+ resl(0) = F3d(0,0);
+ resl(1) = F3d(0,1);
+ resl(2) = F3d(0,2);
+ resl(3) = F3d(1,0);
+ resl(4) = F3d(1,1);
+ resl(5) = F3d(1,2);
+ resl(6) = F3d(2,0);
+ resl(7) = F3d(2,1);
+ resl(8) = F3d(2,2);
+ break;
+
+ case OutputType::outGrp_strain:
+ case OutputType::outGrp_C:
+ case OutputType::outGrp_I1: {
+ // In-plane Cauchy-Green deformation tensor
+ auto C = mat_mul(transpose(F), F);
+
+ // In-plane Green-Lagrange strain tensor
+ auto Eg = 0.50 * (C - Im);
+
+ if (outGrp == OutputType::outGrp_strain) {
+ // resl is used to remap Eg
+ resl(0) = Eg(0,0);
+ resl(1) = Eg(1,1);
+ resl(2) = Eg(2,2);
+ resl(3) = Eg(0,1);
+ resl(4) = Eg(1,2);
+ resl(5) = Eg(2,0);
+
+ } else if (outGrp == OutputType::outGrp_C) {
+ // resl is used to remap C
+ resl(0) = C(0,0);
+ resl(1) = C(1,1);
+ resl(2) = C(2,2);
+ resl(3) = C(0,1);
+ resl(4) = C(1,2);
+ resl(5) = C(2,0);
+
+ } else if (outGrp == OutputType::outGrp_I1) {
+ resl(0) = mat_trace(C, 3);
+ sE(e) = sE(e) + w*resl(0);
+ }
+ } break;
+
+ case OutputType::outGrp_stress: {
+ //dmsg << "outGrp: " << " outGrp_stress";
+ Array S(3,3);
+
+ // 2nd Piola-Kirchhoff stress
+ S(0,0) = Sm(0,0);
+ S(1,1) = Sm(1,0);
+ S(0,1) = Sm(2,0);
+ S(1,0) = S(0,1);
+
+ // Normalizing stress by thickness
+ S = S / ht;
+
+ // 2nd Piola-Kirchhoff stress tensor
+ resl(0) = S(0,0);
+ resl(1) = S(1,1);
+ resl(2) = S(2,2);
+ resl(3) = S(0,1);
+ resl(4) = S(1,2);
+ resl(5) = S(2,0);
+ //dmsg << "resl: " << resl;
+ } break;
+ }
+
+ for (int a = 0; a < lM.eNoN; a++) {
+ int Ac = lM.IEN(a,e);
+ sA(Ac)= sA(Ac) + w*N(a);
+ for (int i = 0; i < m; i++) {
+ sF(i,Ac) = sF(i,Ac) + w*N(a)*resl(i);
+ }
+ }
+ }
+
+ if (!is_zero(Je)) {
+ sE(e) = sE(e) / Je;
+ }
+ }
+
+ resE = sE;
+
+ // Exchange data at the shared nodes across processes
+ all_fun::commu(com_mod, sF);
+ all_fun::commu(com_mod, sA);
+
+ for (int a = 0; a < lM.nNo; a++) {
+ int Ac = lM.gN(a);
+ if (!is_zero(sA(Ac))) {
+ for (int i = 0; i < m; i++) {
+ res(i,a) = res(i,a) + sF(i,Ac) / sA(Ac);
+ }
+ }
+ }
+}
+
//-------
// tpost
//-------
diff --git a/Code/Source/svFSI/post.h b/Code/Source/svFSI/post.h
index 79061175..2fec5150 100644
--- a/Code/Source/svFSI/post.h
+++ b/Code/Source/svFSI/post.h
@@ -25,6 +25,9 @@ void post(Simulation* simulation, const mshType& lM, Array& res, const A
void ppbin2vtk(Simulation* simulation);
+void shl_post(Simulation* simulation, const mshType& lM, const int m, Array& res,
+ Vector& resE, const Array& lD, const int iEq, consts::OutputType outGrp);
+
void tpost(Simulation* simulation, const mshType& lM, const int m, Array& res, Vector& resE, const Array& lD,
const Array& lY, const int iEq, consts::OutputType outGrp);
diff --git a/Code/Source/svFSI/read_files.cpp b/Code/Source/svFSI/read_files.cpp
index 704be82c..fc8d9487 100644
--- a/Code/Source/svFSI/read_files.cpp
+++ b/Code/Source/svFSI/read_files.cpp
@@ -422,8 +422,8 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq,
// Read BCs for shells with triangular elements. Not necessary for
// NURBS elements
//
- if (bc_params->shell_bc_type.defined()) {
- auto ctmp = bc_params->shell_bc_type.value();
+ if (bc_params->cst_shell_bc_type.defined()) {
+ auto ctmp = bc_params->cst_shell_bc_type.value();
if (std::set{"Fixed", "fixed", "Clamped", "clamped"}.count(ctmp)) {
lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_fix));
if (!utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Dir))) {
@@ -1166,8 +1166,12 @@ void read_domain(Simulation* simulation, EquationParameters* eq_params, eqType&
read_cep_domain(simulation, eq_params, domain_params, lEq.dmn[iDmn]);
}
- // Set parameters for a solid material model.
- if ((lEq.dmn[iDmn].phys == EquationType::phys_struct) || (lEq.dmn[iDmn].phys == EquationType::phys_ustruct)) {
+ // Read material/constitutive model parameters for nonlinear
+ // elastodynamics simulations (both solids and shells)
+ //
+ if ( (lEq.dmn[iDmn].phys == EquationType::phys_shell) ||
+ (lEq.dmn[iDmn].phys == EquationType::phys_struct) ||
+ (lEq.dmn[iDmn].phys == EquationType::phys_ustruct)) {
read_mat_model(simulation, eq_params, domain_params, lEq.dmn[iDmn]);
if (utils::is_zero(lEq.dmn[iDmn].stM.Kpen) && lEq.dmn[iDmn].phys == EquationType::phys_struct) {
//err = "Incompressible struct is not allowed. Use "// "penalty method or ustruct"
@@ -1440,7 +1444,16 @@ void read_files(Simulation* simulation, const std::string& file_name)
auto& com_mod = simulation->get_com_mod();
+ #define n_debug_read_files
+ #ifdef debug_read_files
+ DebugMsg dmsg(__func__, com_mod.cm.idcm());
+ dmsg.banner();
+ #endif
+
// Read the solver XML file.
+ #ifdef debug_read_files
+ dmsg << "Read the solver XML file " << " ... ";
+ #endif
if (!com_mod.resetSim) {
simulation->read_parameters(std::string(file_name));
}
@@ -1479,6 +1492,9 @@ void read_files(Simulation* simulation, const std::string& file_name)
simulation->set_module_parameters();
// Read mesh and BCs data.
+ #ifdef debug_read_files
+ dmsg << "Read mesh and BCs data " << " ... ";
+ #endif
read_msh_ns::read_msh(simulation);
// Reading immersed boundary mesh data.
@@ -1502,6 +1518,10 @@ void read_files(Simulation* simulation, const std::string& file_name)
com_mod.nEq = nEq;
com_mod.eq.resize(nEq);
std::for_each(com_mod.eq.begin(),com_mod.eq.end(),[&](eqType& eq){eq.roInf=simulation->roInf;});
+ #ifdef debug_read_files
+ dmsg << "Read equations " << " ... ";
+ dmsg << "nEq: " << nEq;
+ #endif
for (int iEq = 0; iEq < nEq; iEq++) {
auto& eq = com_mod.eq[iEq];
@@ -1573,6 +1593,9 @@ void read_files(Simulation* simulation, const std::string& file_name)
}
}
}
+ #ifdef debug_read_files
+ dmsg << "Done Read equations " << " ";
+ #endif
auto& cep_mod = simulation->get_cep_mod();
@@ -1635,6 +1658,10 @@ void read_files(Simulation* simulation, const std::string& file_name)
com_mod.cplBC.nX = 0;
com_mod.cplBC.xo.clear();
}
+
+ #ifdef debug_read_files
+ dmsg << "Done" << " ";
+ #endif
}
//--------------------------------
@@ -1933,19 +1960,23 @@ void read_mat_model(Simulation* simulation, EquationParameters* eq_params, Domai
}
// If no constitutive model was given use a NeoHookean model.
+ //
+ ConstitutiveModelType cmodel_type;
+ std::string cmodel_str;
+
if (!domain_params->constitutive_model.defined()) {
lDmn.stM.isoType = ConstitutiveModelType::stIso_nHook;
- lDmn.stM.C10 = mu * 0.5;
- return;
- }
+ cmodel_type = ConstitutiveModelType::stIso_nHook;
+ cmodel_str = "neoHookean";
// Get the constitutive model type.
- ConstitutiveModelType cmodel_type;
- auto cmodel_str = domain_params->constitutive_model.type.value();
- try {
- cmodel_type = constitutive_model_name_to_type.at(cmodel_str);
- } catch (const std::out_of_range& exception) {
- throw std::runtime_error("Unknown constitutive model type '" + cmodel_str + ".");
+ } else {
+ cmodel_str = domain_params->constitutive_model.type.value();
+ try {
+ cmodel_type = constitutive_model_name_to_type.at(cmodel_str);
+ } catch (const std::out_of_range& exception) {
+ throw std::runtime_error("Unknown constitutive model type '" + cmodel_str + ".");
+ }
}
// Set material properties for the domain 'lDmn'.
@@ -1972,6 +2003,21 @@ void read_mat_model(Simulation* simulation, EquationParameters* eq_params, Domai
}
}
+ // Check for shell model
+ //
+ // ST91 is the default and the only dilational penalty model for
+ // compressible shell elements. This is set to avoid any square-
+ // root evaulations of the Jacobian during Newton iterations for
+ // satisfying plane-stress condition.
+ //
+ if (lDmn.phys == EquationType::phys_shell) {
+ lDmn.stM.Kpen = kap;
+ if (!incompFlag) {
+ lDmn.stM.volType = ConstitutiveModelType::stVol_ST91;
+ }
+ return;
+ }
+
// Look for dilational penalty model. HGO uses quadratic penalty model.
if (domain_params->dilational_penalty_model.defined()) {
auto model_str = domain_params->dilational_penalty_model.value();
diff --git a/Code/Source/svFSI/read_msh.cpp b/Code/Source/svFSI/read_msh.cpp
index 8194f4ac..982d220a 100644
--- a/Code/Source/svFSI/read_msh.cpp
+++ b/Code/Source/svFSI/read_msh.cpp
@@ -1237,6 +1237,7 @@ void read_msh(Simulation* simulation)
if (mesh.eType != consts::ElementType::NRB && mesh.eType != consts::ElementType::TRI3) {
throw std::runtime_error("Shell elements must be triangles or C1-NURBS.");
}
+
if (mesh.eType == consts::ElementType::NRB) {
for (int i = 0; i < com_mod.nsd-1; i++) {
if (mesh.bs[i].p <= 1) {
@@ -1244,6 +1245,13 @@ void read_msh(Simulation* simulation)
}
}
}
+
+ if (mesh.eType == consts::ElementType::TRI3) {
+ if (!com_mod.cm.seq()) {
+ throw std::runtime_error("Shells with linear triangles should be run sequentially.");
+ }
+ }
+
}
}
@@ -1675,33 +1683,41 @@ void read_msh(Simulation* simulation)
// Read contact model parameters.
//
- // [NOTE] Implement this, no tests yet.
- //
if (!com_mod.resetSim) {
com_mod.iCntct = false;
- //lPM => list%get(ctmp, "Contact model")
- /*
- IF (ASSOCIATED(lPM)) THEN
- iCntct = .TRUE.
- SELECT CASE (TRIM(ctmp))
- CASE ("penalty")
- cntctM%cType = cntctM_penalty
- lPtr => lPM%get(cntctM%k,
- 2 "Penalty constant (k)", 1, ll=0._RKIND)
- lPtr => lPM%get(cntctM%h,
- 2 "Desired separation (h)", 1, lb=0._RKIND)
- lPtr => lPM%get(cntctM%c,
- 2 "Closest gap to activate penalty (c)", 1, lb=0._RKIND)
- IF (cntctM%c .LT. cntctM%h) err =
- 2 "Choose c > h for proper contact penalization"
- lPtr => lPM%get(cntctM%al,
- 2 "Min norm of face normals (alpha)",1,lb=0._RKIND,
- 3 ub=1._RKIND)
- CASE DEFAULT
- err = "Undefined contact model"
- END SELECT
- END IF
- */
+ auto& cntctM = com_mod.cntctM;
+ auto& contact_params = simulation->parameters.contact_parameters;
+
+ if (contact_params.model.defined()) {
+ auto contact_model = contact_params.model.value();
+ com_mod.iCntct = true;
+
+ try {
+ cntctM.cType = consts::contact_model_name_to_type.at(contact_model);
+ } catch (const std::out_of_range& exception) {
+ throw std::runtime_error("Unknown contact model '" + contact_model + "'.");
+ }
+
+ switch (cntctM.cType) {
+
+ case consts::ContactModelType::cntctM_penalty:
+ cntctM.k = contact_params.penalty_constant.value();
+ cntctM.h = contact_params.desired_separation.value();
+ cntctM.c = contact_params.closest_gap_to_activate_penalty.value();
+ cntctM.al = contact_params.min_norm_of_face_normals.value();
+
+ if (cntctM.c < cntctM.h) {
+ throw std::runtime_error("The contact Closest_gap_to_activate_penalty " + std::to_string(cntctM.c) +
+ " must be > the desired separation " + std::to_string(cntctM.h) + ".");
+ }
+ break;
+
+ default:
+ throw std::runtime_error("Contact model '" + contact_model + "' is not implemented.");
+ break;
+ }
+
+ }
}
}
diff --git a/Code/Source/svFSI/set_equation_props.h b/Code/Source/svFSI/set_equation_props.h
index 84cd9183..17f8dbfb 100644
--- a/Code/Source/svFSI/set_equation_props.h
+++ b/Code/Source/svFSI/set_equation_props.h
@@ -455,6 +455,7 @@ SetEquationPropertiesMapType set_equation_props = {
using namespace consts;
auto& com_mod = simulation->get_com_mod();
lEq.phys = consts::EquationType::phys_shell;
+ com_mod.shlEq = true;
propL[0][0] = PhysicalProperyType::solid_density;
propL[1][0] = PhysicalProperyType::damping;
@@ -467,8 +468,18 @@ SetEquationPropertiesMapType set_equation_props = {
read_domain(simulation, eq_params, lEq, propL);
- nDOP = {3,1,1,0};
- outPuts = {OutputType::out_displacement, OutputType::out_velocity, OutputType::out_integ};
+ nDOP = {9,1,0,0};
+ outPuts = {
+ OutputType::out_displacement,
+ OutputType::out_stress,
+ OutputType::out_strain,
+ OutputType::out_jacobian,
+ OutputType::out_defGrad,
+ OutputType::out_velocity,
+ OutputType::out_integ,
+ OutputType::out_CGstrain,
+ OutputType::out_CGInv1
+ };
// Set solver parameters.
read_ls(simulation, eq_params, SolverType::lSolver_CG, lEq);
diff --git a/Code/Source/svFSI/set_material_props.h b/Code/Source/svFSI/set_material_props.h
index 96ebde15..61bc74f8 100644
--- a/Code/Source/svFSI/set_material_props.h
+++ b/Code/Source/svFSI/set_material_props.h
@@ -125,4 +125,23 @@ SeMaterialPropertiesMapType set_material_props = {
lDmn.stM.bfs = params.bfs.value();
} },
+//---------------------------//
+// stIso_LS //
+//---------------------------//
+// Lee-Sacks material model.
+//
+{consts::ConstitutiveModelType::stIso_LS, [](DomainParameters* domain_params, double mu, double kap, double lam,
+ dmnType& lDmn) -> void
+{
+ lDmn.stM.isoType = consts::ConstitutiveModelType::stIso_LS;
+ auto& params = domain_params->constitutive_model.lee_sacks;
+
+ lDmn.stM.a = params.a.value();
+ lDmn.stM.a0 = params.a0.value();
+ lDmn.stM.b1 = params.b1.value();
+ lDmn.stM.b2 = params.b2.value();
+ lDmn.stM.mu0 = params.mu0.value();
+
+} },
+
};
diff --git a/Code/Source/svFSI/set_output_props.h b/Code/Source/svFSI/set_output_props.h
index 6cd1fc6f..1272c47f 100644
--- a/Code/Source/svFSI/set_output_props.h
+++ b/Code/Source/svFSI/set_output_props.h
@@ -16,6 +16,7 @@ using OutputProps = std::tuple;
//------------------
// output_props_map
//------------------
+// Reproduces Fortran READOUTPUTS.
//
std::map output_props_map =
{
@@ -25,12 +26,19 @@ std::map output_props_map =
{OutputType::out_absVelocity, std::make_tuple(OutputType::outGrp_absV, 0, nsd, "Absolute_velocity") },
{OutputType::out_acceleration, std::make_tuple(OutputType::outGrp_A, 0, nsd, "Acceleration") },
{OutputType::out_cauchy, std::make_tuple(OutputType::outGrp_cauchy, 0, com_mod.nsymd, "Cauchy_stress") },
+
+ {OutputType::out_CGInv1, std::make_tuple(OutputType::out_CGInv1, 0, 1, "CG_Strain_Trace") },
+ {OutputType::out_CGstrain, std::make_tuple(OutputType::outGrp_C, 0, com_mod.nsymd, "CG_Strain") },
+
{OutputType::out_defGrad, std::make_tuple(OutputType::outGrp_F, 0, nsd*nsd, "Def_grad") },
{OutputType::out_displacement, std::make_tuple(OutputType::outGrp_D, 0, nsd, "Displacement") },
- {OutputType::out_divergence, std::make_tuple(OutputType::outGrp_divV, 0, 1, "Divergence") },
+ {OutputType::out_divergence, std::make_tuple(OutputType::outGrp_divV, 0, 1, "Divergence") },
{OutputType::out_energyFlux, std::make_tuple(OutputType::outGrp_eFlx, 0, nsd, "Energy_flux") },
- {OutputType::out_fibAlign, std::make_tuple(OutputType::outGrp_fA, 0, 1, "Fiber_alignment") },
+
+ {OutputType::out_fibAlign, std::make_tuple(OutputType::outGrp_fA, 0, 1, "Fiber_alignment") },
{OutputType::out_fibDir, std::make_tuple(OutputType::outGrp_fN, 0, nsd, "Fiber_direction") },
+ {OutputType::out_fibStrn, std::make_tuple(OutputType::outGrp_fS, 0, 1, "Fiber_shortening") },
+
{OutputType::out_heatFlux, std::make_tuple(OutputType::outGrp_hFlx, 0, nsd, "Heat_flux") },
{OutputType::out_integ, std::make_tuple(OutputType::outGrp_I, 0, 1, nsd == 2 ? "Area" : "Volume") },
{OutputType::out_jacobian, std::make_tuple(OutputType::outGrp_J, 0, 1, "Jacobian") },
diff --git a/Code/Source/svFSI/shells.cpp b/Code/Source/svFSI/shells.cpp
index d896310b..0f995f4e 100644
--- a/Code/Source/svFSI/shells.cpp
+++ b/Code/Source/svFSI/shells.cpp
@@ -4,16 +4,1671 @@
#include "shells.h"
+#include "all_fun.h"
+#include "consts.h"
+#include "lhsa.h"
+#include "mat_fun.h"
+#include "mat_models.h"
#include "nn.h"
+#include "utils.h"
+
+#ifdef WITH_TRILINOS
+#include "trilinos_linear_solver.h"
+#endif
namespace shells {
+//-----------------
+// construct_shell
+//-----------------
+// This routines is for solving nonlinear shell mechanics problem
+// using finite elements and IGA.
+//
+// Reproduces Fortran CONSTRUCT_SHELL
+//
+void construct_shell(ComMod& com_mod, const mshType& lM, const Array& Ag,
+ const Array& Yg, const Array& Dg)
+{
+ using namespace consts;
+
+ #define n_debug_construct_shell
+ #ifdef debug_construct_shell
+ DebugMsg dmsg(__func__, com_mod.cm.idcm());
+ dmsg.banner();
+ dmsg << "lM.nFn: " << lM.nFn;
+ dmsg << "lM.nFs: " << lM.nFs;
+ dmsg << "lM.eNoN: " << lM.eNoN;
+ #endif
+
+ const int nsd = com_mod.nsd;
+ const int tDof = com_mod.tDof;
+ const int dof = com_mod.dof;
+ const int cEq = com_mod.cEq;
+ const auto& eq = com_mod.eq[cEq];
+ auto cDmn = com_mod.cDmn;
+
+ int eNoN = lM.eNoN;
+
+ if (lM.eType == ElementType::TRI3) {
+ eNoN = 2 * eNoN;
+ }
+
+ int nFn = lM.nFn;
+ if (nFn == 0) {
+ nFn = 1;
+ }
+
+ // Initialize tensor operations
+ mat_fun::ten_init(2);
+
+ // SHELLS: dof = nsd
+ //
+ Vector ptr(eNoN);
+ Array xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN),
+ bfl(nsd,eNoN), fN(3,nFn), lR(dof,eNoN);
+ Array3 lK(dof*dof,eNoN,eNoN);
+
+ // Loop over all elements of mesh
+ //
+ for (int e = 0; e < lM.nEl; e++) {
+ // Update domain and proceed if domain phys and eqn phys match
+ cDmn = all_fun::domain(com_mod, lM, cEq, e);
+ auto cPhys = eq.dmn[cDmn].phys;
+ if (cPhys != EquationType::phys_shell) {
+ continue;
+ }
+
+#if 0
+ dmsg << " " << " ";
+ dmsg << " " << " ";
+ dmsg << "-------------------------------------" << " ";
+ dmsg << "--------------- e: " << e+1;
+ dmsg << "-------------------------------------" << " ";
+#endif
+
+ // Create local copies
+ xl = 0.0;
+ al = 0.0;
+ yl = 0.0;
+ dl = 0.0;
+ bfl = 0.0;
+
+ for (int a = 0; a < eNoN; a++) {
+ //dmsg << "----- a: " << a;
+ int Ac = -1;
+
+ if (a < lM.eNoN) {
+ Ac = lM.IEN(a,e);
+ ptr(a) = Ac;
+ } else {
+ int b = a - lM.eNoN;
+ //dmsg << "b: " << b;
+ Ac = lM.eIEN(b,e);
+ ptr(a) = Ac;
+ //dmsg << "Ac: " << Ac;
+ if (Ac == -1) {
+ continue;
+ }
+ }
+
+ //dmsg << "Ac: " << Ac;
+
+ for (int i = 0; i < nsd; i++) {
+ xl(i,a) = com_mod.x(i,Ac);
+ bfl(i,a) = com_mod.Bf(i,Ac);
+ }
+
+ for (int i = 0; i < tDof; i++) {
+ al(i,a) = Ag(i,Ac);
+ dl(i,a) = Dg(i,Ac);
+ yl(i,a) = Yg(i,Ac);
+ }
+ }
+
+ //dmsg << "lM.fN.size(): " << lM.fN.size();
+ if (lM.fN.size() != 0) {
+ for (int iFn = 0; iFn < nFn; iFn++) {
+ for (int i = 0; i < 3; i++) {
+ fN(i,iFn) = lM.fN(i+nsd*iFn,e);
+ }
+ }
+ } else {
+ fN = 0.0;
+ }
+
+ #ifdef debug_construct_shell
+ dmsg << "-------------------" << "-------------------";
+ dmsg << "ptr: " << ptr;
+ dmsg << "-------------------" << "-------------------";
+ #endif
+
+ // Constant strain triangles, no numerical integration
+ //
+ if (lM.eType == ElementType::TRI3) {
+ shell_cst(com_mod, lM, e, eNoN, nFn, fN, al, yl, dl, xl, bfl, ptr);
+
+ } else {
+ lR = 0.0;
+ lK = 0.0;
+
+ // Update shape functions for NURBS elements
+ //if (lM.eType .EQ. eType_NRB) CALL NRBNNX(lM, e)
+
+ // Gauss integration
+ for (int g = 0; g < lM.nG; g++) {
+ shell_3d(com_mod, lM, g, eNoN, nFn, fN, al, yl, dl, xl, bfl, lR, lK);
+ }
+ }
+
+ //if (e+1 == 19) {
+ //exit(0);
+ //}
+
+ // Assembly
+#ifdef WITH_TRILINOS
+ if (eq.assmTLS) {
+ trilinos_doassem_(const_cast(eNoN), const_cast(ptr.data()), lK.data(), lR.data());
+ } else {
+#endif
+ lhsa_ns::do_assem(com_mod, eNoN, ptr, lK, lR);
+#ifdef WITH_TRILINOS
+ }
+#endif
+
+ } // e: loop
+
+}
+
+//----------
+// shell_3d
+//----------
+// Construct shell mechanics for higher order elements/NURBS
+//
+void shell_3d(ComMod& com_mod, const mshType& lM, const int g, const int eNoN,
+ const int nFn, const Array& fN,
+ const Array& al, const Array& yl, const Array& dl, const Array& xl,
+ const Array& bfl, Array& lR, Array3& lK)
+{
+ std::cout << "========== shell_3d ==========" << std::endl;
+ std::cout << "[shell_3d] g: " << g << std::endl;
+
+ using namespace consts;
+ using namespace mat_fun;
+
+ const int nsd = com_mod.nsd;
+ const int dof = com_mod.dof;
+ int cEq = com_mod.cEq;
+ auto& eq = com_mod.eq[cEq];
+ auto cDmn = com_mod.cDmn;
+ auto& dmn = eq.dmn[cDmn];
+ const double dt = com_mod.dt;
+
+ // Define parameters
+ double rho = eq.dmn[cDmn].prop.at(PhysicalProperyType::solid_density);
+ double dmp = dmn.prop.at(PhysicalProperyType::damping);
+ double ht = eq.dmn[cDmn].prop.at(PhysicalProperyType::shell_thickness);
+ Vector fb({dmn.prop.at(PhysicalProperyType::f_x), dmn.prop.at(PhysicalProperyType::f_y),
+ dmn.prop.at(PhysicalProperyType::f_z)});
+ double amd = eq.am * rho + eq.af * eq.gam * dt * dmp;
+ double afl = eq.af * eq.beta * dt * dt;
+
+ int i = eq.s;
+ int j = i + 1;
+ int k = j + 1;
+
+ // Get the reference configuration
+ auto x0 = xl;
+
+ // Get the current configuration
+ //
+ Array xc(3,eNoN);
+
+ for (int a = 0; a < eNoN; a++) {
+ xc(0,a) = x0(0,a) + dl(i,a);
+ xc(1,a) = x0(1,a) + dl(j,a);
+ xc(2,a) = x0(2,a) + dl(k,a);
+ }
+
+ // Define shape functions and their derivatives at Gauss point
+ //
+ Vector N;
+ Array Nx, Nxx;
+
+ /* [TODO] Nurbs are not supported.
+ if (lM.eType == ElementType::eType_NRB) {
+ N = lM.N.rcol(g);
+ Nx = lM.Nx.rslice(g);
+ Nxx = lM.Nxx.rslice(g);
+ } else {
+ N = lM.fs(0).N.rcol(g);
+ Nx = lM.fs(0).Nx.rslice(g);
+ Nxx = lM.fs(0).Nxx.rslice(g);
+ }
+ */
+ N = lM.fs[0].N.rcol(g);
+ Nx = lM.fs[0].Nx.rslice(g);
+ Nxx = lM.fs[0].Nxx.rslice(g);
+
+//=====================================================================
+// TODO: Might have to call GNNxx for Jacobian transformation. Check
+// formulation again.
+//
+// CALL GNNxx(2, eNoN, 2, lM.fs(0).Nx(:,:,g), lM.fs(0).Nxx(:,:,g),
+// xl, Nx, Nxx)
+//
+//=====================================================================
+
+ // Compute preliminaries on the reference configuration
+ // Covariant and contravariant bases (reference config)
+ //
+ Array aCov0(3,2), aCnv0(3,2);
+ Vector nV0(3);
+ nn::gnns(nsd, lM.eNoN, Nx, x0, nV0, aCov0, aCnv0);
+ double Jac0 = sqrt(utils::norm(nV0));
+ nV0 = nV0 / Jac0;
+
+ // Second derivatives for computing curvature coeffs. (ref. config)
+ //
+ Array3 r0_xx(2,2,3);
+
+ for (int a = 0; a < eNoN; a++) {
+ for (int i = 0; i < 3; i++) {
+ r0_xx(0,0,i) = r0_xx(0,0,i) + Nxx(0,a)*x0(i,a);
+ r0_xx(1,1,i) = r0_xx(1,1,i) + Nxx(1,a)*x0(i,a);
+ r0_xx(0,1,i) = r0_xx(0,1,i) + Nxx(2,a)*x0(i,a);
+ }
+ }
+
+ for (int i = 0; i < 3; i++) {
+ r0_xx(1,0,i) = r0_xx(0,1,i);
+ }
+
+ // Compute metric tensor and curvature coefficients (ref. config)
+ //
+ double aa_0[2][2]{}, bb_0[2][2]{};
+
+ for (int l = 0; l < nsd; l++) {
+ aa_0[0][0] = aa_0[0][0] + aCov0(l,0)*aCov0(l,0);
+ aa_0[0][1] = aa_0[0][1] + aCov0(l,0)*aCov0(l,1);
+ aa_0[1][0] = aa_0[1][0] + aCov0(l,1)*aCov0(l,0);
+ aa_0[1][1] = aa_0[1][1] + aCov0(l,1)*aCov0(l,1);
+
+ bb_0[0][0] = bb_0[0][0] + r0_xx(0,0,l)*nV0(l);
+ bb_0[0][1] = bb_0[0][1] + r0_xx(0,1,l)*nV0(l);
+ bb_0[1][0] = bb_0[1][0] + r0_xx(1,0,l)*nV0(l);
+ bb_0[1][1] = bb_0[1][1] + r0_xx(1,1,l)*nV0(l);
+ }
+
+ // Compute fiber orientation in curvature coordinates
+ //
+ Array fNa0(2,nFn);
+
+ for (int iFn = 0; iFn < nFn; iFn++) {
+ for (int l = 0; l < 3; l++) {
+ fNa0(0,iFn) = fNa0(0,iFn) + fN(l,iFn)*aCnv0(l,0);
+ fNa0(1,iFn) = fNa0(1,iFn) + fN(l,iFn)*aCnv0(l,1);
+ }
+ }
+
+ // Now compute preliminaries on the current configuration
+ // Covariant and contravariant bases (current/spatial config)
+ //
+ Array aCov(3,2), aCnv(3,2);
+ Vector nV(3);
+ nn::gnns(nsd, eNoN, Nx, xc, nV, aCov, aCnv);
+ double Jac = sqrt(utils::norm(nV));
+ nV = nV / Jac;
+
+ // Second derivatives for computing curvature coeffs. (cur. config)
+ Array3 r_xx(2,2,3);
+
+ for (int a = 0; a < eNoN; a++) {
+ r_xx(0,0,i) = r_xx(0,0,i) + Nxx(0,a)*xc(i,a);
+ r_xx(1,1,i) = r_xx(1,1,i) + Nxx(1,a)*xc(i,a);
+ r_xx(0,1,i) = r_xx(0,1,i) + Nxx(2,a)*xc(i,a);
+ }
+
+ for (int i = 0; i < 3; i++) {
+ r_xx(1,0,i) = r_xx(0,1,i);
+ }
+
+ // Compute metric tensor and curvature coefficients (cur. config)
+ //
+ // [TODO] are these dimensions correct? is nsd = 3?
+ //
+ double aa_x[2][2]{}, bb_x[2][2]{};
+
+ for (int l = 0; l < nsd; l++) {
+ aa_x[0][0] = aa_x[0][0] + aCov(l,0)*aCov(l,0);
+ aa_x[0][1] = aa_x[0][1] + aCov(l,0)*aCov(l,1);
+ aa_x[1][0] = aa_x[1][0] + aCov(l,1)*aCov(l,0);
+ aa_x[1][1] = aa_x[1][1] + aCov(l,1)*aCov(l,1);
+
+ bb_x[0][0] = bb_x[0][0] + r_xx(0,0,l)*nV(l);
+ bb_x[0][1] = bb_x[0][1] + r_xx(0,1,l)*nV(l);
+ bb_x[1][0] = bb_x[1][0] + r_xx(1,0,l)*nV(l);
+ bb_x[1][1] = bb_x[1][1] + r_xx(1,1,l)*nV(l);
+ }
+
+ // Compute stress resultants by integrating 2nd Piola Kirchhoff
+ // stress and elasticity tensors through the shell thickness. These
+ // resultants are computed in Voigt notation.
+ //
+ Array3 Dm(3,3,3);
+ Array Sm(3,2);
+ double lam3;
+
+ shl_strs_res(com_mod, dmn, nFn, fNa0, aa_0, aa_x, bb_0, bb_x, lam3, Sm, Dm);
+
+ // Variation in the membrane strain
+ //
+ Array3 Bm(3,3,eNoN);
+
+ for (int a = 0; a < eNoN; a++) {
+ Bm(0,0,a) = Nx(0,a)*aCov(0,0);
+ Bm(0,1,a) = Nx(0,a)*aCov(1,0);
+ Bm(0,2,a) = Nx(0,a)*aCov(2,0);
+
+ Bm(1,0,a) = Nx(1,a)*aCov(0,1);
+ Bm(1,1,a) = Nx(1,a)*aCov(1,1);
+ Bm(1,2,a) = Nx(1,a)*aCov(2,1);
+
+ Bm(2,0,a) = Nx(1,a)*aCov(0,0) + Nx(0,a)*aCov(0,1);
+ Bm(2,1,a) = Nx(1,a)*aCov(1,0) + Nx(0,a)*aCov(1,1);
+ Bm(2,2,a) = Nx(1,a)*aCov(2,0) + Nx(0,a)*aCov(2,1);
+ }
+
+ // Variation in the bending strain
+ // dB = -(B1 + B2) du; B1 = N_xx * n;
+ // B2 = (r_xx Nm M1 Nx - r_xx N M2 Nx)
+ //
+ // Second derivatives of the position vector (current)
+ //
+ Array Kc(3,3);
+
+ for (int i = 0; i < 3; i++) {
+ Kc(0,i) = r_xx(0,0,i);
+ Kc(1,i) = r_xx(1,1,i);
+ Kc(2,i) = r_xx(0,1,i) + r_xx(1,0,i);
+ }
+
+ // N matrix
+ auto Nm = mat_id(3) - mat_dyad_prod(nV, nV, 3);
+ Nm = Nm / Jac;
+
+ // M1, M2 matrices
+ //
+ Array Mm(3,3);
+ Array3 KNmMm(3,3,2);
+
+ for (int l = 0; l < 2; l++) {
+ Mm = 0.0;
+ Mm(0,1) = -aCov(2,l);
+ Mm(0,2) = aCov(1,l);
+ Mm(1,2) = -aCov(0,l);
+
+ // Skew-symmetric
+ Mm(1,0) = -Mm(0,1);
+ Mm(2,0) = -Mm(0,2);
+ Mm(2,1) = -Mm(1,2);
+
+ KNmMm.set_slice(l, mat_mul(Kc, mat_mul(Nm, Mm)));
+ }
+
+ // Define variation in bending strain tensor (Bb), Voigt notation
+ //
+ Array3 Bb(3,3,eNoN);
+
+ for (int a = 0; a < eNoN; a++) {
+ for (int i = 0; i < 3; i++) {
+ Bb(0,i,a) = -Nxx(0,a)*nV(i);
+ Bb(1,i,a) = -Nxx(1,a)*nV(i);
+ Bb(2,i,a) = -Nxx(2,a)*nV(i)*2.0;
+ }
+ }
+
+ for (int a = 0; a < eNoN; a++) {
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ Bb(i,j,a) = Bb(i,j,a) + Nx(0,a)*KNmMm(i,j,1) - Nx(1,a)*KNmMm(i,j,0);
+ }
+ }
+ }
+
+ // Contribution to tangent matrices: Dm * Bm, Dm*Bb
+ //
+ Array3 D0Bm(3,3,eNoN), D1Bm(3,3,eNoN), D1Bb(3,3,eNoN), D2Bb(3,3,eNoN);
+
+ for (int a = 0; a < eNoN; a++) {
+ D0Bm.set_slice(a, mat_mul(Dm.rslice(0), Bm.rslice(a)));
+ D1Bm.set_slice(a, mat_mul(Dm.rslice(1), Bm.rslice(a)));
+ D1Bb.set_slice(a, mat_mul(Dm.rslice(1), Bb.rslice(a)));
+ D2Bb.set_slice(a, mat_mul(Dm.rslice(2), Bb.rslice(a)));
+ }
+
+ // Acceleration and mass damping at the integration point
+ //
+ auto ud = -fb;
+
+ for (int a = 0; a < eNoN; a++) {
+ ud(0) = ud(0) + N(a)*(rho*(al(i,a)-bfl(0,a)) + dmp*yl(i,a));
+ ud(1) = ud(1) + N(a)*(rho*(al(j,a)-bfl(1,a)) + dmp*yl(j,a));
+ ud(2) = ud(2) + N(a)*(rho*(al(k,a)-bfl(2,a)) + dmp*yl(k,a));
+ }
+
+ // Local residue
+ //
+ auto w = lM.w(g)*Jac0;
+ auto wh = w*ht;
+
+ for (int a = 0; a < eNoN; a++) {
+ double BmS = Bm(0,0,a)*Sm(0,0) + Bm(1,0,a)*Sm(1,0) + Bm(2,0,a)*Sm(2,0);
+ double BbS = Bb(0,0,a)*Sm(0,1) + Bb(1,0,a)*Sm(1,1) + Bb(2,0,a)*Sm(2,1);
+ lR(0,a) = lR(0,a) + wh*N(a)*ud(0) + w*(BmS + BbS);
+
+ BmS = Bm(0,1,a)*Sm(0,0) + Bm(1,1,a)*Sm(1,0) + Bm(2,1,a)*Sm(2,0);
+ BbS = Bb(0,1,a)*Sm(0,1) + Bb(1,1,a)*Sm(1,1) + Bb(2,1,a)*Sm(2,1);
+ lR(1,a) = lR(1,a) + wh*N(a)*ud(1) + w*(BmS + BbS);
+
+ BmS = Bm(0,2,a)*Sm(0,0) + Bm(1,2,a)*Sm(1,0) + Bm(2,2,a)*Sm(2,0);
+ BbS = Bb(0,2,a)*Sm(0,1) + Bb(1,2,a)*Sm(1,1) + Bb(2,2,a)*Sm(2,1);
+ lR(2,a) = lR(2,a) + wh*N(a)*ud(2) + w*(BmS + BbS);
+ }
+
+ // Local stiffness
+ //
+ amd = wh*amd;
+ afl = w*afl;
+
+ for (int b = 0; b < eNoN; b++) {
+ for (int a = 0; a < eNoN; a++) {
+ // Contribution from inertia and geometric stiffness
+ double NxSNx = Nx(0,a)*Nx(0,b)*Sm(0,0) + Nx(1,a)*Nx(1,b)*Sm(1,0) + Nx(0,a)*Nx(1,b)*Sm(2,0) +
+ Nx(1,a)*Nx(0,b)*Sm(2,0);
+ double T1 = amd*N(a)*N(b) + afl*NxSNx;
+
+ lK(0,a,b) = lK(0,a,b) + T1;
+ lK(dof+1,a,b) = lK(dof+1,a,b) + T1;
+ lK(2*dof+2,a,b) = lK(2*dof+2,a,b) + T1;
+
+ // Contribution from material stiffness
+ double BmDBm = Bm(0,0,a)*D0Bm(0,0,b) + Bm(1,0,a)*D0Bm(1,0,b) + Bm(2,0,a)*D0Bm(2,0,b);
+ double BmDBb = Bm(0,0,a)*D1Bb(0,0,b) + Bm(1,0,a)*D1Bb(1,0,b) + Bm(2,0,a)*D1Bb(2,0,b);
+ double BbDBm = Bb(0,0,a)*D1Bm(0,0,b) + Bb(1,0,a)*D1Bm(1,0,b) + Bb(2,0,a)*D1Bm(2,0,b);
+ double BbDBb = Bb(0,0,a)*D2Bb(0,0,b) + Bb(1,0,a)*D2Bb(1,0,b) + Bb(2,0,a)*D2Bb(2,0,b);
+ lK(0,a,b) = lK(0,a,b) + afl*(BmDBm + BmDBb + BbDBm + BbDBb);
+
+ BmDBm = Bm(0,0,a)*D0Bm(0,1,b) + Bm(1,0,a)*D0Bm(1,1,b) + Bm(2,0,a)*D0Bm(2,1,b);
+ BmDBb = Bm(0,0,a)*D1Bb(0,1,b) + Bm(1,0,a)*D1Bb(1,1,b) + Bm(2,0,a)*D1Bb(2,1,b);
+ BbDBm = Bb(0,0,a)*D1Bm(0,1,b) + Bb(1,0,a)*D1Bm(1,1,b) + Bb(2,0,a)*D1Bm(2,1,b);
+ BbDBb = Bb(0,0,a)*D2Bb(0,1,b) + Bb(1,0,a)*D2Bb(1,1,b) + Bb(2,0,a)*D2Bb(2,1,b);
+ lK(1,a,b) = lK(1,a,b) + afl*(BmDBm + BmDBb + BbDBm + BbDBb);
+
+ BmDBm = Bm(0,0,a)*D0Bm(0,2,b) + Bm(1,0,a)*D0Bm(1,2,b) + Bm(2,0,a)*D0Bm(2,2,b);
+ BmDBb = Bm(0,0,a)*D1Bb(0,2,b) + Bm(1,0,a)*D1Bb(1,2,b) + Bm(2,0,a)*D1Bb(2,2,b);
+ BbDBm = Bb(0,0,a)*D1Bm(0,2,b) + Bb(1,0,a)*D1Bm(1,2,b) + Bb(2,0,a)*D1Bm(2,2,b);
+ BbDBb = Bb(0,0,a)*D2Bb(0,2,b) + Bb(1,0,a)*D2Bb(1,2,b) + Bb(2,0,a)*D2Bb(2,2,b);
+ lK(2,a,b) = lK(2,a,b) + afl*(BmDBm + BmDBb + BbDBm + BbDBb);
+
+ BmDBm = Bm(0,1,a)*D0Bm(0,0,b) + Bm(1,1,a)*D0Bm(1,0,b) + Bm(2,1,a)*D0Bm(2,0,b);
+ BmDBb = Bm(0,1,a)*D1Bb(0,0,b) + Bm(1,1,a)*D1Bb(1,0,b) + Bm(2,1,a)*D1Bb(2,0,b);
+ BbDBm = Bb(0,1,a)*D1Bm(0,0,b) + Bb(1,1,a)*D1Bm(1,0,b) + Bb(2,1,a)*D1Bm(2,0,b);
+ BbDBb = Bb(0,1,a)*D2Bb(0,0,b) + Bb(1,1,a)*D2Bb(1,0,b) + Bb(2,1,a)*D2Bb(2,0,b);
+ lK(dof+0,a,b) = lK(dof+0,a,b) + afl*(BmDBm + BmDBb + BbDBm + BbDBb);
+
+ BmDBm = Bm(0,1,a)*D0Bm(0,1,b) + Bm(1,1,a)*D0Bm(1,1,b) + Bm(2,1,a)*D0Bm(2,1,b);
+ BbDBm = Bb(0,1,a)*D1Bm(0,1,b) + Bb(1,1,a)*D1Bm(1,1,b) + Bb(2,1,a)*D1Bm(2,1,b);
+ BbDBb = Bb(0,1,a)*D2Bb(0,1,b) + Bb(1,1,a)*D2Bb(1,1,b) + Bb(2,1,a)*D2Bb(2,1,b);
+ lK(dof+1,a,b) = lK(dof+1,a,b) + afl*(BmDBm + BmDBb + BbDBm + BbDBb);
+
+ BmDBm = Bm(0,1,a)*D0Bm(0,2,b) + Bm(1,1,a)*D0Bm(1,2,b) + Bm(2,1,a)*D0Bm(2,2,b);
+ BmDBb = Bm(0,1,a)*D1Bb(0,2,b) + Bm(1,1,a)*D1Bb(1,2,b) + Bm(2,1,a)*D1Bb(2,2,b);
+ BbDBm = Bb(0,1,a)*D1Bm(0,2,b) + Bb(1,1,a)*D1Bm(1,2,b) + Bb(2,1,a)*D1Bm(2,2,b);
+ BbDBb = Bb(0,1,a)*D2Bb(0,2,b) + Bb(1,1,a)*D2Bb(1,2,b) + Bb(2,1,a)*D2Bb(2,2,b);
+ lK(dof+2,a,b) = lK(dof+2,a,b) + afl*(BmDBm + BmDBb + BbDBm + BbDBb);
+
+ BmDBm = Bm(0,2,a)*D0Bm(0,0,b) + Bm(1,2,a)*D0Bm(1,0,b) + Bm(2,2,a)*D0Bm(2,0,b);
+ BmDBb = Bm(0,2,a)*D1Bb(0,0,b) + Bm(1,2,a)*D1Bb(1,0,b) + Bm(2,2,a)*D1Bb(2,0,b);
+ BbDBm = Bb(0,2,a)*D1Bm(0,0,b) + Bb(1,2,a)*D1Bm(1,0,b) + Bb(2,2,a)*D1Bm(2,0,b);
+ BbDBb = Bb(0,2,a)*D2Bb(0,0,b) + Bb(1,2,a)*D2Bb(1,0,b) + Bb(2,2,a)*D2Bb(2,0,b);
+ lK(2*dof+0,a,b) = lK(2*dof+0,a,b) + afl*(BmDBm + BmDBb + BbDBm + BbDBb);
+
+ BmDBm = Bm(0,2,a)*D0Bm(0,1,b) + Bm(1,2,a)*D0Bm(1,1,b) + Bm(2,2,a)*D0Bm(2,1,b);
+ BmDBb = Bm(0,2,a)*D1Bb(0,1,b) + Bm(1,2,a)*D1Bb(1,1,b) + Bm(2,2,a)*D1Bb(2,1,b);
+ BbDBm = Bb(0,2,a)*D1Bm(0,1,b) + Bb(1,2,a)*D1Bm(1,1,b) + Bb(2,2,a)*D1Bm(2,1,b);
+ BbDBb = Bb(0,2,a)*D2Bb(0,1,b) + Bb(1,2,a)*D2Bb(1,1,b) + Bb(2,2,a)*D2Bb(2,1,b);
+ lK(2*dof+1,a,b) = lK(2*dof+1,a,b) + afl*(BmDBm + BmDBb + BbDBm + BbDBb);
+
+ BmDBm = Bm(0,2,a)*D0Bm(0,2,b) + Bm(1,2,a)*D0Bm(1,2,b) + Bm(2,2,a)*D0Bm(2,2,b);
+ BmDBb = Bm(0,2,a)*D1Bb(0,2,b) + Bm(1,2,a)*D1Bb(1,2,b) + Bm(2,2,a)*D1Bb(2,2,b);
+ BbDBm = Bb(0,2,a)*D1Bm(0,2,b) + Bb(1,2,a)*D1Bm(1,2,b) + Bb(2,2,a)*D1Bm(2,2,b);
+ BbDBb = Bb(0,2,a)*D2Bb(0,2,b) + Bb(1,2,a)*D2Bb(1,2,b) + Bb(2,2,a)*D2Bb(2,2,b);
+ lK(2*dof+2,a,b) = lK(2*dof+2,a,b) + afl*(BmDBm + BmDBb + BbDBm + BbDBb);
+ }
+ }
+}
+
+//----------------
+// shell_bend_cst
+//----------------
+// This subroutine computes bending strain, Eb, and its variational
+// derivative, Bb, for CST elements
+//
+// Reproduces Fortran SHELLBENDCST.
+//
+void shell_bend_cst(ComMod& com_mod, const mshType& lM, const int e, const Vector& ptr,
+ Array& x0, Array& xc, double bb_0[2][2], double bb_x[2][2],
+ Array3& Bb, const bool vflag)
+{
+ using namespace consts;
+ using namespace mat_fun;
+ using namespace utils;
+
+ #define n_debug_shell_bend_cst
+ #ifdef debug_shell_bend_cst
+ DebugMsg dmsg(__func__, com_mod.cm.idcm());
+ dmsg.banner();
+ dmsg << "e: " << e;
+ #endif
+
+ int cEq = com_mod.cEq;
+ auto& eq = com_mod.eq[cEq];
+ auto cDmn = com_mod.cDmn;
+ auto& dmn = eq.dmn[cDmn];
+
+ // Define parameters
+ double rho = eq.dmn[cDmn].prop.at(PhysicalProperyType::solid_density);
+ double dmp = dmn.prop.at(PhysicalProperyType::damping);
+
+ int nsd = com_mod.nsd;
+ int eNoN = 2 * lM.eNoN;
+
+ // Boundary element check
+ bool bFlag = false;
+
+ for (int j = lM.eNoN; j < eNoN; j++) {
+ if (ptr(j) == -1) {
+ bFlag = true;
+ break;
+ }
+ }
+
+ // Edge vectors of the main element (reference config)
+ //
+ Array a0(3,6);
+
+ for (int i = 0; i < 3; i++) {
+ a0(i,0) = x0(i,2) - x0(i,1);
+ a0(i,1) = x0(i,0) - x0(i,2);
+ a0(i,2) = x0(i,1) - x0(i,0);
+ }
+
+ // Edge vectors of the main element (current config)
+ //
+ Array a(3,6);
+
+ for (int i = 0; i < 3; i++) {
+ a(i,0) = xc(i,2) - xc(i,1);
+ a(i,1) = xc(i,0) - xc(i,2);
+ a(i,2) = xc(i,1) - xc(i,0);
+ }
+
+ // Covariant and contravariant bases in reference config
+ Array tmpA(3,3);
+
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ tmpA(i,j) = x0(i,j);
+ }
+ }
+
+ Array aCov0(3,2), aCnv0(3,2);
+ Vector nV0(3);
+ nn::gnns(nsd, lM.eNoN, lM.Nx.rslice(0), tmpA, nV0, aCov0, aCnv0);
+ double Jac0 = sqrt(utils::norm(nV0));
+ nV0 = nV0 / Jac0;
+
+ // Covariant and contravariant bases in current config
+ //
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ tmpA(i,j) = xc(i,j);
+ }
+ }
+
+ Array aCov(3,2), aCnv(3,2);
+ Vector nV(3);
+ nn::gnns(nsd, lM.eNoN, lM.Nx.rslice(0), tmpA, nV, aCov, aCnv);
+ double Jac = sqrt(norm(nV));
+ nV = nV / Jac;
+
+ // Update the position vector of the `artificial' or `ghost' nodes
+ // depending on the boundary condition.
+ //
+ if (bFlag) {
+ for (int j = lM.eNoN; j < eNoN; j++) {
+ if (ptr(j) != -1) {
+ continue;
+ }
+ int i = j - lM.eNoN;
+ int p = i - 1;
+
+ if (i == 0) {
+ p = 2;
+ }
+
+ // Reference config
+ // eI = eI0 = aI0/|aI0| (reference config)
+ //
+ auto aIi = 1.0 / sqrt(norm(a0.col(i)));
+ auto eI = a0.col(i) * aIi;
+
+ // nI = nI0 = eI0 x n0 (reference config)
+ //
+ Vector nI(3);
+ nI(0) = eI(1)*nV0(2) - eI(2)*nV0(1);
+ nI(1) = eI(2)*nV0(0) - eI(0)*nV0(2);
+ nI(2) = eI(0)*nV0(1) - eI(1)*nV0(0);
+
+ // xJ = xI + 2(nI \ctimes nI)aP
+ //
+ auto nInI = mat_dyad_prod(nI, nI, 3);
+ x0(0,j) = 2.0 * (nInI(0,0)*a0(0,p) + nInI(0,1)*a0(1,p) + nInI(0,2)*a0(2,p)) + x0(0,i);
+ x0(1,j) = 2.0 * (nInI(1,0)*a0(0,p) + nInI(1,1)*a0(1,p) + nInI(1,2)*a0(2,p)) + x0(1,i);
+ x0(2,j) = 2.0 * (nInI(2,0)*a0(0,p) + nInI(2,1)*a0(1,p) + nInI(2,2)*a0(2,p)) + x0(2,i);
+
+ // Current config
+ // eI = aI/|aI| (current config)
+ //
+ aIi = 1.0 / sqrt(utils::norm(a.col(i)));
+ eI = a.col(i)*aIi;
+
+ // nI = eI x n (currnt config)
+ //
+ nI(0) = eI(1)*nV(2) - eI(2)*nV(1);
+ nI(1) = eI(2)*nV(0) - eI(0)*nV(2);
+ nI(2) = eI(0)*nV(1) - eI(1)*nV(0);
+
+ // xJ = xI + 2(nI \ctimes nI)aP
+ //
+ nInI = mat_dyad_prod(nI, nI, 3);
+ xc(0,j) = 2.0*(nInI(0,0)*a(0,p) + nInI(0,1)*a(1,p) + nInI(0,2)*a(2,p)) + xc(0,i);
+ xc(1,j) = 2.0*(nInI(1,0)*a(0,p) + nInI(1,1)*a(1,p) + nInI(1,2)*a(2,p)) + xc(1,i);
+ xc(2,j) = 2.0*(nInI(2,0)*a(0,p) + nInI(2,1)*a(1,p) + nInI(2,2)*a(2,p)) + xc(2,i);
+
+ if (utils::btest(lM.sbc(i,e),enum_int(BoundaryConditionType::bType_fix))) {
+ for (int i = 0; i < 3; i++) {
+ xc(i,j) = x0(i,j);
+ }
+ }
+ }
+ }
+
+ // Edge vector of surrounding nodes (reference config)
+ //
+ for (int i = 0; i < 3; i++) {
+ a0(i,3) = x0(i,3) - x0(i,1);
+ a0(i,4) = x0(i,4) - x0(i,2);
+ a0(i,5) = x0(i,5) - x0(i,0);
+ }
+
+ // Edge vector of surrounding nodes (current config)
+ //
+ for (int i = 0; i < 3; i++) {
+ a(i,3) = xc(i,3) - xc(i,1);
+ a(i,4) = xc(i,4) - xc(i,2);
+ a(i,5) = xc(i,5) - xc(i,0);
+ }
+
+ // a.gCnv (reference config)
+ //
+ Array adg0(3,3);
+ adg0(0,0) = norm(a0.col(3), aCnv0.col(0)); // xi_4
+ adg0(0,1) = norm(a0.col(4), aCnv0.col(0)); // xi_5
+ adg0(0,2) = norm(a0.col(5), aCnv0.col(0)); // xi_6
+
+ adg0(1,0) = norm(a0.col(3), aCnv0.col(1)); // eta_4
+ adg0(1,1) = norm(a0.col(4), aCnv0.col(1)); // eta_5
+ adg0(1,2) = norm(a0.col(5), aCnv0.col(1)); // eta_6
+
+ adg0(2,0) = norm(a0.col(3), nV0); // z_4
+ adg0(2,1) = norm(a0.col(4), nV0); // z_5
+ adg0(2,2) = norm(a0.col(5), nV0); // z_6
+
+ // a.gCnv (current config)
+ //
+ Array adg(3,3);
+ adg(0,0) = norm(a.col(3), aCnv.col(0)); // xi_4
+ adg(0,1) = norm(a.col(4), aCnv.col(0)); // xi_5
+ adg(0,2) = norm(a.col(5), aCnv.col(0)); // xi_6
+
+ adg(1,0) = norm(a.col(3), aCnv.col(1)); // eta_4
+ adg(1,1) = norm(a.col(4), aCnv.col(1)); // eta_5
+ adg(1,2) = norm(a.col(5), aCnv.col(1)); // eta_6
+
+ adg(2,0) = norm(a.col(3), nV); // z_4
+ adg(2,1) = norm(a.col(4), nV); // z_5
+ adg(2,2) = norm(a.col(5), nV); // z_6
+
+ // Xi matrix (reference config)
+ //
+ auto xi0 = adg0;
+ xi0(0,2) = adg0(0,2) + 1.0; // xi_6
+ xi0(1,0) = adg0(1,0) + 1.0; // eta_4
+
+ // Xi matrix (current config)
+ //
+ auto xi = adg;
+ xi(0,2) = adg(0,2) + 1.0; // xi_6
+ xi(1,0) = adg(1,0) + 1.0; // eta_4
+
+ // Tmat and inverse (reference config)
+ //
+ Array Tm0(3,3);
+
+ for (int i = 0; i < 3; i++) {
+ Tm0(i,0) = xi0(0,i)*(xi0(0,i) - 1.0); // xi**2 - xi
+ Tm0(i,1) = xi0(1,i)*(xi0(1,i) - 1.0); // eta**2 - eta
+ Tm0(i,2) = xi0(0,i)*xi0(1,i); // xi * eta
+ }
+
+ Tm0 = mat_inv(Tm0, 3);
+
+ // Tmat and inverse (current config)
+ //
+ Array Tm(3,3);
+
+ for (int i = 0; i < 3; i++) {
+ Tm(i,0) = xi(0,i)*(xi(0,i) - 1.0); // xi**2 - xi
+ Tm(i,1) = xi(1,i)*(xi(1,i) - 1.0); // eta**2 - eta
+ Tm(i,2) = xi(0,i)*xi(1,i); // xi * eta
+ }
+
+ Tm = mat_inv(Tm, 3);
+
+ // v = Inv(T) * z (reference config)
+ //
+ Vector v0(3);
+
+ v0(0) = Tm0(0,0)*xi0(2,0) + Tm0(0,1)*xi0(2,1) + Tm0(0,2)*xi0(2,2);
+ v0(1) = Tm0(1,0)*xi0(2,0) + Tm0(1,1)*xi0(2,1) + Tm0(1,2)*xi0(2,2);
+ v0(2) = Tm0(2,0)*xi0(2,0) + Tm0(2,1)*xi0(2,1) + Tm0(2,2)*xi0(2,2);
+
+ // Curvature coefficients (ref. config)
+ bb_0[0][0] = 2.0 * v0(0);
+ bb_0[1][1] = 2.0 * v0(1);
+ bb_0[0][1] = v0(2);
+ bb_0[1][0] = bb_0[0][1];
+
+ // v = Inv(T) * z (current config)
+ //
+ double v[3];
+ v[0] = Tm(0,0)*xi(2,0) + Tm(0,1)*xi(2,1) + Tm(0,2)*xi(2,2);
+ v[1] = Tm(1,0)*xi(2,0) + Tm(1,1)*xi(2,1) + Tm(1,2)*xi(2,2);
+ v[2] = Tm(2,0)*xi(2,0) + Tm(2,1)*xi(2,1) + Tm(2,2)*xi(2,2);
+
+ // Curvature coefficients (current config)
+ //
+ bb_x[0][0] = 2.0*v[0];
+ bb_x[1][1] = 2.0*v[1];
+ bb_x[0][1] = v[2];
+ bb_x[1][0] = bb_x[0][1];
+
+ if (!vflag) {
+ Bb = 0.0;
+ return;
+ }
+
+ // Now compute variation in bending strain
+ //
+ // B1 bar
+ //
+ Array B1b(3,6);
+
+ for (int i = 0; i < 3; i++) {
+ B1b(i,0) = -Tm(i,0) * ((2.0*xi(0,0)-1.0)*v[0] + xi(1,0)*v[2]);
+ B1b(i,1) = -Tm(i,0) * ((2.0*xi(1,0)-1.0)*v[1] + xi(0,0)*v[2]);
+
+ B1b(i,2) = -Tm(i,1) * ((2.0*xi(0,1)-1.0)*v[0] + xi(1,1)*v[2]);
+ B1b(i,3) = -Tm(i,1) * ((2.0*xi(1,1)-1.0)*v[1] + xi(0,1)*v[2]);
+
+ B1b(i,4) = -Tm(i,2) * ((2.0*xi(0,2)-1.0)*v[0] + xi(1,2)*v[2]);
+ B1b(i,5) = -Tm(i,2) * ((2.0*xi(1,2)-1.0)*v[1] + xi(0,2)*v[2]);
+ }
+
+ // H1
+ //
+ Array H1(6,18);
+
+ for (int i = 0; i < 3; i++) {
+ H1(0, i) = aCnv(i,0)*adg(1,0);
+ H1(1, i) = aCnv(i,1)*adg(1,0);
+ H1(2, i) = aCnv(i,0)*adg(1,1);
+ H1(3, i) = aCnv(i,1)*adg(1,1);
+ H1(4, i) = aCnv(i,0)*adg(1,2);
+ H1(5, i) = aCnv(i,1)*adg(1,2);
+
+ H1(0, i+3 ) = -aCnv(i,0)*adg(0,0);
+ H1(1, i+3 ) = -aCnv(i,1)*adg(0,0);
+ H1(2, i+3 ) = -aCnv(i,0)*adg(0,1);
+ H1(3, i+3 ) = -aCnv(i,1)*adg(0,1);
+ H1(4, i+3 ) = -aCnv(i,0)*adg(0,2);
+ H1(5, i+3 ) = -aCnv(i,1)*adg(0,2);
+
+ H1(0,i+9) = aCnv(i,0);
+ H1(1,i+9) = aCnv(i,1);
+
+ H1(2,i+12) = aCnv(i,0);
+ H1(3,i+12) = aCnv(i,1);
+
+ H1(4,i+15) = aCnv(i,0);
+ H1(5,i+15) = aCnv(i,1);
+ }
+
+ // H2
+ Array H2(18,18);
+ H2( 0, 3) = -1.0;
+ H2( 1, 4) = -1.0;
+ H2( 2, 5) = -1.0;
+ H2( 3, 6) = -1.0;
+ H2( 4, 7) = -1.0;
+ H2( 5, 8) = -1.0;
+ H2( 6, 0) = -1.0;
+ H2( 7, 1) = -1.0;
+ H2( 8, 2) = -1.0;
+
+ H2( 0, 6) = 1.0;
+ H2( 1, 7) = 1.0;
+ H2( 2, 8) = 1.0;
+ H2( 3, 0) = 1.0;
+ H2( 4, 1) = 1.0;
+ H2( 5, 2) = 1.0;
+ H2( 6, 3) = 1.0;
+ H2( 7, 4) = 1.0;
+ H2( 8, 5) = 1.0;
+
+ H2(9, 3) = -1.0;
+ H2(10, 4) = -1.0;
+ H2(11, 5) = -1.0;
+ H2(12, 6) = -1.0;
+ H2(13, 7) = -1.0;
+ H2(14, 8) = -1.0;
+ H2(15, 0) = -1.0;
+ H2(16, 1) = -1.0;
+ H2(17, 2) = -1.0;
+
+ for (int i = 9; i < 18; i++) {
+ H2(i,i) = 1.0;
+ }
+
+ // N matrix
+ auto Nm = mat_id(3) - mat_dyad_prod(nV, nV, 3);
+
+ // M1, M2 matrices
+ //
+ Array3 Mm(3,3,2);
+
+ for (int i = 0; i < 2; i++) {
+ Mm(0,1,i) = -aCov(2,i);
+ Mm(0,2,i) = aCov(1,i);
+ Mm(1,2,i) = -aCov(0,i);
+
+ // Skew-symmetric
+ Mm(1,0,i) = -Mm(0,1,i);
+ Mm(2,0,i) = -Mm(0,2,i);
+ Mm(2,1,i) = -Mm(1,2,i);
+ }
+
+ // H3 matrix
+ //
+ Array H3(3,18);
+ tmpA = mat_mul(Nm, Mm.rslice(0));
+ tmpA = -tmpA / Jac;
+
+ H3(0,0) = a(0,3)*tmpA(0,0) + a(1,3)*tmpA(1,0) + a(2,3)*tmpA(2,0);
+ H3(0,1) = a(0,3)*tmpA(0,1) + a(1,3)*tmpA(1,1) + a(2,3)*tmpA(2,1);
+ H3(0,2) = a(0,3)*tmpA(0,2) + a(1,3)*tmpA(1,2) + a(2,3)*tmpA(2,2);
+
+ H3(1,0) = a(0,4)*tmpA(0,0) + a(1,4)*tmpA(1,0) + a(2,4)*tmpA(2,0);
+ H3(1,1) = a(0,4)*tmpA(0,1) + a(1,4)*tmpA(1,1) + a(2,4)*tmpA(2,1);
+ H3(1,2) = a(0,4)*tmpA(0,2) + a(1,4)*tmpA(1,2) + a(2,4)*tmpA(2,2);
+
+ H3(2,0) = a(0,5)*tmpA(0,0) + a(1,5)*tmpA(1,0) + a(2,5)*tmpA(2,0);
+ H3(2,1) = a(0,5)*tmpA(0,1) + a(1,5)*tmpA(1,1) + a(2,5)*tmpA(2,1);
+ H3(2,2) = a(0,5)*tmpA(0,2) + a(1,5)*tmpA(1,2) + a(2,5)*tmpA(2,2);
+
+ tmpA = mat_mul(Nm, Mm.rslice(1));
+ tmpA = -tmpA / Jac;
+
+ H3(0,3) = a(0,3)*tmpA(0,0) + a(1,3)*tmpA(1,0) + a(2,3)*tmpA(2,0);
+ H3(0,4) = a(0,3)*tmpA(0,1) + a(1,3)*tmpA(1,1) + a(2,3)*tmpA(2,1);
+ H3(0,5) = a(0,3)*tmpA(0,2) + a(1,3)*tmpA(1,2) + a(2,3)*tmpA(2,2);
+
+ H3(1,3) = a(0,4)*tmpA(0,0) + a(1,4)*tmpA(1,0) + a(2,4)*tmpA(2,0);
+ H3(1,4) = a(0,4)*tmpA(0,1) + a(1,4)*tmpA(1,1) + a(2,4)*tmpA(2,1);
+ H3(1,5) = a(0,4)*tmpA(0,2) + a(1,4)*tmpA(1,2) + a(2,4)*tmpA(2,2);
+
+ H3(2,3) = a(0,5)*tmpA(0,0) + a(1,5)*tmpA(1,0) + a(2,5)*tmpA(2,0);
+ H3(2,4) = a(0,5)*tmpA(0,1) + a(1,5)*tmpA(1,1) + a(2,5)*tmpA(2,1);
+ H3(2,5) = a(0,5)*tmpA(0,2) + a(1,5)*tmpA(1,2) + a(2,5)*tmpA(2,2);
+
+ for (int i = 0; i < 3; i++) {
+ H3(0,i+9) = nV(i);
+ H3(1,i+12) = nV(i);
+ H3(2,i+15) = nV(i);
+ }
+
+ // Variation in bending strain (Bb = -2*(B1b*H1*H2 + Tinv*H3*H2))
+ auto H1H2 = mat_mul(H1, H2);
+ auto Bb1 = mat_mul(B1b, H1H2);
+ auto H3H2 = mat_mul(H3, H2);
+ Bb1 = -2.0 * (Bb1 + mat_mul(Tm, H3H2));
+ Bb.set_values(Bb1);
+ //Bb = RESHAPE(Bb1, SHAPE(Bb))
+
+ // Update Bb for boundary elements
+ if (bFlag) {
+ std::vector lFix = {false, false, false};
+ auto Im = mat_id(3);
+
+ for (int j = lM.eNoN; j < eNoN; j++) {
+ if (ptr(j) != -1) {
+ continue;
+ }
+
+ int i = j - lM.eNoN;
+ int p = i - 1;
+ int f = i + 1;
+
+ if (i == 0) {
+ p = 2;
+ }
+
+ if (i == 2) {
+ f = 0;
+ }
+
+ double aIi, cI;
+ Vector eI(3), nI(3);
+ Array nInI(3,3), eIeI(3,3), eIaP(3,3);
+
+ if (utils::btest(lM.sbc(i,e),enum_int(BoundaryConditionType::bType_fix))) {
+ // eI = eI0 = aI0/|aI0| (reference config)
+ aIi = 1.0 / sqrt(norm(a0.col(i)));
+ eI = a0.col(i) * aIi;
+
+ // nI = nI0 = eI0 x n0 (reference config)
+ nI(0) = eI(1)*nV0(2) - eI(2)*nV0(1);
+ nI(1) = eI(2)*nV0(0) - eI(0)*nV0(2);
+ nI(2) = eI(0)*nV0(1) - eI(1)*nV0(0);
+ nInI = mat_dyad_prod(nI, nI, 3);
+
+ } else {
+ // eI = aI/|aI| (current config)
+ aIi = 1.0 / sqrt(norm(a.col(i)));
+ eI = a.col(i)*aIi;
+
+ // nI = eI x n (currnt config)
+ //
+ Vector nI(3);
+
+ nI(0) = eI(1)*nV(2) - eI(2)*nV(1);
+ nI(1) = eI(2)*nV(0) - eI(0)*nV(2);
+ nI(2) = eI(0)*nV(1) - eI(1)*nV(0);
+
+ cI = norm(a.col(i),a.col(p))*aIi*aIi;
+ nInI = mat_dyad_prod(nI, nI, 3);
+ eIeI = mat_dyad_prod(eI, eI, 3);
+ eIaP = mat_dyad_prod(eI, a.col(p), 3);
+ }
+
+ // Update Bb now
+ // Free boundary conditions: assumed that the `artificial'
+ // triangle is always located in the plane of the main element
+ // of the patch
+ //
+ if (utils::btest(lM.sbc(i,e),enum_int(BoundaryConditionType::bType_free))) {
+
+ // E_I
+ //
+ if (!lFix[i]) {
+ tmpA = -Im + 2.0 * eIeI;
+ Bb.rslice(i) = Bb.rslice(i) + mat_mul(Bb.rslice(j), tmpA);
+ }
+
+ // E_P
+ //
+ if (!lFix[p]) {
+ tmpA = -2.0 *(cI*Im - 2.0 *cI*eIeI + aIi*eIaP);
+ Bb.rslice(p) = Bb.rslice(p) + mat_mul(Bb.rslice(j), tmpA);
+ }
+
+ // E_F
+ //
+ if (!lFix[f]) {
+ tmpA = 2.0 * ((1.0-cI)*Im - (1.0-2.0*cI)*eIeI -aIi*eIaP);
+ Bb.rslice(f) = Bb.rslice(f) + mat_mul(Bb.rslice(j), tmpA);
+ }
+
+ Bb.rslice(j) = 0.0;
+
+ // Hinged boundary conditions: a special case of simple support
+ // in which no translation displacements are allowed.
+
+ } else if (utils::btest(lM.sbc(i,e),enum_int(BoundaryConditionType::bType_hing))) {
+
+ // E_I
+ if (!lFix[i]) {
+ tmpA = -Im + 2.0*eIeI;
+ Bb.rslice(i) = Bb.rslice(i) + mat_mul(Bb.rslice(j), tmpA);
+ }
+
+ lFix[p] = true;
+ lFix[f] = true;
+ Bb.rslice(p) = 0.0;
+ Bb.rslice(f) = 0.0;
+ Bb.rslice(j) = 0.0;
+
+ // Fixed boundary condition: no displacements and no rotations
+ // are allowed.
+ //
+ } else if (utils::btest(lM.sbc(i,e),enum_int(BoundaryConditionType::bType_fix))) {
+
+ if (!lFix[i]) {
+ tmpA = Im - 2.0*nInI;
+ Bb.rslice(i) = Bb.rslice(i) + mat_mul(Bb.rslice(j), tmpA);
+ }
+
+ lFix[p] = true;
+ lFix[f] = true;
+ Bb.rslice(f) = 0.0;
+ Bb.rslice(p) = 0.0;
+
+ // Symmetric BCs (need to be verified)
+ //
+ } else if (utils::btest(lM.sbc(i,e),enum_int(BoundaryConditionType::bType_symm))) {
+ if (!lFix[i]) {
+ tmpA = Im - 2.0*nInI;
+ Bb.rslice(i) = Bb.rslice(i) + mat_mul(Bb.rslice(j), tmpA);
+ }
+
+ tmpA.rcol(0) = eI;
+ tmpA.rcol(1) = nV;
+ tmpA.rcol(2) = nI;
+
+ Bb.rslice(j) = 0.0;
+ Bb.rslice(f) = mat_mul(Bb.rslice(f), tmpA);
+ Bb.rslice(p) = mat_mul(Bb.rslice(p), tmpA);
+
+ for (int i = 0; i < 3; i++) {
+ Bb(i,2,f) = 0.0;
+ Bb(i,2,p) = 0.0;
+ }
+ }
+ }
+ }
+}
+
+//----------
+// shell_bf
+//----------
+// Set follower pressure load/net traction on shells. The traction
+// on shells is treated as body force and the subroutine is called
+// from BF.f
+//
+// Reproduces Fortran SHELLBF.
+//
+void shell_bf(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx,
+ const Array& dl, const Array& xl, const Array& tfl, Array& lR, Array3& lK)
+{
+ using namespace consts;
+
+ const int nsd = com_mod.nsd;
+ const int dof = com_mod.dof;
+ int cEq = com_mod.cEq;
+ auto& eq = com_mod.eq[cEq];
+ auto cDmn = com_mod.cDmn;
+ auto& dmn = eq.dmn[cDmn];
+ const double dt = com_mod.dt;
+
+ double afl = eq.af * eq.beta * dt * dt;
+
+ int i = eq.s;
+ int j = i + 1;
+ int k = j + 1;
+
+ // Get the current configuration and traction vector
+ //
+ double tfn = 0.0;
+ Array xc(3,eNoN);
+ // [NOTE] This is a hack for enabling 'tfl' to be used as a vector in the Fortran.
+ auto tfl_data = tfl.data();
+
+ for (int a = 0; a < eNoN; a++) {
+ xc(0,a) = xl(0,a) + dl(i,a);
+ xc(1,a) = xl(1,a) + dl(j,a);
+ xc(2,a) = xl(2,a) + dl(k,a);
+
+ tfn = tfn + N(a)*tfl_data[a];
+ }
+
+ double wl = w * tfn;
+
+ // Covariant and contravariant bases in current config
+ Array gCov(3,2), gCnv(3,2);
+ Vector nV(3);
+ nn::gnns(nsd, eNoN, Nx, xc, nV, gCov, gCnv);
+
+ // Local residue
+ for (int a = 0; a < eNoN; a++) {
+ lR(0,a) = lR(0,a) - wl*N(a)*nV(0);
+ lR(1,a) = lR(1,a) - wl*N(a)*nV(1);
+ lR(2,a) = lR(2,a) - wl*N(a)*nV(2);
+ }
+
+ // Local stiffness: mass matrix and stiffness contribution due to
+ // follower traction load
+ //
+ double T1 = afl*wl*0.50;
+
+ for (int b = 0; b < eNoN; b++) {
+ for (int a = 0; a < eNoN; a++) {
+ auto lKp = gCov.rcol(0)*(N(b)*Nx(1,a) - N(a)*Nx(1,b)) -
+ gCov.rcol(1)*(N(b)*Nx(0,a) - N(a)*Nx(0,b));
+
+ lK(1,a,b) = lK(1,a,b) - T1*lKp(2);
+ lK(2,a,b) = lK(2,a,b) + T1*lKp(1);
+
+ lK(dof+0,a,b) = lK(dof+0,a,b) + T1*lKp(2);
+ lK(dof+2,a,b) = lK(dof+2,a,b) - T1*lKp(0);
+
+ lK(2*dof+0,a,b) = lK(2*dof+0,a,b) - T1*lKp(1);
+ lK(2*dof+1,a,b) = lK(2*dof+1,a,b) + T1*lKp(0);
+ }
+ }
+}
+
+//-----------
+// shell_cst
+//-----------
+// Construct shell mechanics for constant strain triangle elements
+//
+// Note that for triangular elements, eNoN=6 and lM.eNoN=3
+//
+// Reproduces Fortran SHELLCST.
+//
+void shell_cst(ComMod& com_mod, const mshType& lM, const int e, const int eNoN, const int nFn, const Array& fN,
+ const Array& al, const Array