From f3b9d9cf6a016c473266131ddca529fdb80490a4 Mon Sep 17 00:00:00 2001 From: Dave Parker Date: Wed, 5 Jul 2023 21:01:30 -0700 Subject: [PATCH] Add new shell implementation 69 (#75) * Add shells and material models. * Add follower forces function. * Add processing Lee-Sacks parameters. * Add Lee-Sacks material model, xml parameters for it, fix some shell processing bugs. * Add debugging messages, remove left-over debugging code in fft, causing a crash. * Fix indexing bugs, set ATOL to be a bit smaller than Fortran. * Post processing for shells. * Add thowing an exception if trying to run in parallel with tri3 shells. * Add shell contact. * Add some debug, check for not Contact parameter. * Add missing include DbgMsg. --- Code/Source/svFSI/Array.h | 17 +- Code/Source/svFSI/Array3.h | 24 + Code/Source/svFSI/ComMod.h | 11 +- Code/Source/svFSI/Parameters.cpp | 140 +- Code/Source/svFSI/Parameters.h | 53 +- Code/Source/svFSI/VtkData.cpp | 16 +- Code/Source/svFSI/baf_ini.cpp | 6 +- Code/Source/svFSI/bf.cpp | 6 +- Code/Source/svFSI/consts.cpp | 8 + Code/Source/svFSI/consts.h | 24 +- Code/Source/svFSI/contact.cpp | 95 +- Code/Source/svFSI/contact.h | 2 +- Code/Source/svFSI/distribute.cpp | 7 +- Code/Source/svFSI/eq_assem.cpp | 3 +- Code/Source/svFSI/fft.cpp | 10 +- Code/Source/svFSI/lhsa.cpp | 8 +- Code/Source/svFSI/main.cpp | 19 +- Code/Source/svFSI/mat_fun.cpp | 59 + Code/Source/svFSI/mat_fun.h | 3 + Code/Source/svFSI/mat_models.cpp | 566 ++++++++ Code/Source/svFSI/mat_models.h | 6 + Code/Source/svFSI/nn.cpp | 6 +- Code/Source/svFSI/post.cpp | 436 ++++++ Code/Source/svFSI/post.h | 3 + Code/Source/svFSI/read_files.cpp | 72 +- Code/Source/svFSI/read_msh.cpp | 66 +- Code/Source/svFSI/set_equation_props.h | 15 +- Code/Source/svFSI/set_material_props.h | 19 + Code/Source/svFSI/set_output_props.h | 12 +- Code/Source/svFSI/shells.cpp | 1775 +++++++++++++++++++++++- Code/Source/svFSI/shells.h | 25 +- Code/Source/svFSI/vtk_xml.cpp | 45 +- Code/Source/svFSILS/lhs.cpp | 1 + 33 files changed, 3427 insertions(+), 131 deletions(-) diff --git a/Code/Source/svFSI/Array.h b/Code/Source/svFSI/Array.h index e3145a2..ed56259 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 8619c26..67816db 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 ee2981e..03f18c9 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 8dbb567..78454a9 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 d5caf8e..737d9ab 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 3221a3a..bf61687 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 59226ed..d720ab6 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 8f16176..259eb19 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 171bb64..a36fc34 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 df0fa06..db62ab2 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 2166f8c..0815d19 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 d0f5fb6..47b6fb7 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 d1ecce0..764e177 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/lhsa.cpp b/Code/Source/svFSI/lhsa.cpp index bf1b5ad..899a1a4 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 accf4fa..90c7c94 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 f8a8db8..bfa1634 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 6531a9b..e168b03 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 c858022..c43a292 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 0b09164..5250f24 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 0959da0..8b1fd56 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 acf1580..f2a4868 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 7906117..2fec515 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 704be82..fc8d948 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 8194f4a..982d220 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 84cd918..17f8dbf 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 96ebde1..61bc74f 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 6cd1fc6..1272c47 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 d896310..1536082 100644 --- a/Code/Source/svFSI/shells.cpp +++ b/Code/Source/svFSI/shells.cpp @@ -4,16 +4,1667 @@ #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" 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_(eNoN, ptr, lK, lR); + } 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& yl, const Array& dl, const Array& xl, + const Array& bfl, const Vector& ptr) +{ + 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; + + #define n_debug_shell_cst + #ifdef debug_shell_cst + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "lM.nFn: " << lM.nFn; + dmsg << "lM.nFs: " << lM.nFs; + dmsg << "dof: " << dof; + dmsg << "eNoN: " << eNoN; + #endif + + // 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; + + #ifdef debug_shell_cst + dmsg << "rho: " << rho; + dmsg << "dmp: " << dmp; + dmsg << "ht: " << ht; + dmsg << "i: " << i; + dmsg << "j: " << j; + dmsg << "k: " << k; + #endif + + // 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); + } + + auto Nx = lM.Nx.slice(0); + + // Covariant and contravariant bases in reference 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); + } + } + + Array aCov0(3,2), aCnv0(3,2); + Vector nV0(3); + + nn::gnns(nsd, lM.eNoN, Nx, tmpX, nV0, aCov0, aCnv0); + //CALL GNNS(lM%eNoN, Nx, tmpX, nV0, aCov0, aCnv0) + + double Jac0 = sqrt(utils::norm(nV0)); + nV0 = nV0 / Jac0; + + // Covariant and contravariant bases in current config + + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < lM.eNoN; j++) { + tmpX(i,j) = xc(i,j); + } + } + + Array aCov(3,2), aCnv(3,2); + Vector nV(3); + + nn::gnns(nsd, lM.eNoN, Nx, tmpX, nV, aCov, aCnv); + + double Jac = sqrt(utils::norm(nV)); + nV = nV / Jac; + + // Compute metric tensor in reference and current config + // + double aa_0[2][2] = {}; + double aa_x[2][2] = {}; + + for (int g = 0; g < nsd; g++) { + aa_0[0][0] = aa_0[0][0] + aCov0(g,0)*aCov0(g,0); + aa_0[0][1] = aa_0[0][1] + aCov0(g,0)*aCov0(g,1); + aa_0[1][0] = aa_0[1][0] + aCov0(g,1)*aCov0(g,0); + aa_0[1][1] = aa_0[1][1] + aCov0(g,1)*aCov0(g,1); + + aa_x[0][0] = aa_x[0][0] + aCov(g,0)*aCov(g,0); + aa_x[0][1] = aa_x[0][1] + aCov(g,0)*aCov(g,1); + aa_x[1][0] = aa_x[1][0] + aCov(g,1)*aCov(g,0); + aa_x[1][1] = aa_x[1][1] + aCov(g,1)*aCov(g,1); + } + + // Compute fiber orientation in curvature coordinates + // + Array fNa0(2,nFn); + //dmsg << "nFn: " << nFn; + //dmsg << "fN: " << fN; + //dmsg << "aCnv0: " << aCnv0; + + 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); + } + } + //dmsg << "fNa0: " << fNa0; + + // Define variation in membrane strain only for the main element + // + Array3 Bm(3,3,lM.eNoN); + + for (int a = 0; a < lM.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); + } + //dmsg << " " << " "; + //dmsg << "Bm: " << Bm; + //exit(0); + + // For the boundary elements, zero-out Bm for fixed/clamped BC. + // + std::vector setIt = {false, false, false}; + int a = lM.eNoN; + + while (a < eNoN) { + //dmsg << "----- a: " << a+1; + if (ptr(a) == -1) { + int b = a - lM.eNoN; + //dmsg << "b: " << b+1; + if (utils::btest(lM.sbc(b,e),enum_int(BoundaryConditionType::bType_fix))) { + //dmsg << "set to true " << b+1; + setIt[b] = true; + } + } + + a = a + 1; + } + + for (int a = 0; a < lM.eNoN; a++) { + if (setIt[a]) { + for (int b = 0; b < lM.eNoN; b++) { + if (a == b) { + continue; + } + Bm.rslice(b) = 0.0; + //dmsg << "Zero slice b: " << b; + } + } + } + //dmsg << " " << " "; + //dmsg << "Bm: " << Bm; + //exit(0); + + // Compute curvature coefficients for bending strain and its + // variation for CST elements + // + double bb_0[2][2], bb_x[2][2]; + Array3 Bb(3,3,eNoN); + shell_bend_cst(com_mod, lM, e, ptr, x0, xc, bb_0, bb_x, Bb, true); + //CALL SHELLBENDCST(lM, e, ptr, x0, xc, bb_0, bb_x, Bb, .TRUE.) + //dmsg << "Bb: " << Bb; + + // Compute stress resultants by integrating 2nd Piola Kirchhoff + // stress and elasticity tensors through the shell thickness. These + // resultants are computed in Voigt notation. + // + /* + dmsg << "aa_0: " << aa_0[0][0]; + dmsg << " : " << aa_0[0][1]; + dmsg << " : " << aa_0[1][0]; + dmsg << " : " << aa_0[1][1]; + dmsg << "aa_x: " << aa_x[0][0]; + dmsg << " : " << aa_x[0][1]; + dmsg << " : " << aa_x[1][0]; + dmsg << " : " << aa_x[1][1]; + */ + 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); + + /* + dmsg << " " << " "; + dmsg << "Bm: " << Bm; + dmsg << "Bb: " << Bb; + + dmsg << " " << " "; + dmsg << "Sm: " << Sm; + dmsg << "Dm: " << Dm; + exit(0); + */ + + // Contribution to tangent matrices: Dm * Bm, Dm*Bb + // + Array3 D0Bm(3,3,lM.eNoN), D1Bm(3,3,lM.eNoN); + + for (int a = 0; a < lM.eNoN; a++) { + D0Bm.rslice(a) = mat_fun::mat_mul(Dm.rslice(0), Bm.rslice(a)); + D1Bm.rslice(a) = mat_fun::mat_mul(Dm.rslice(1), Bm.rslice(a)); + } + + Array3 D1Bb(3,3,eNoN), D2Bb(3,3,eNoN); + + for (int a = 0; a < eNoN; a++) { + D1Bb.rslice(a) = mat_fun::mat_mul(Dm.rslice(1), Bb.rslice(a)); + D2Bb.rslice(a) = mat_fun::mat_mul(Dm.rslice(2), Bb.rslice(a)); + } + + // Contribution to residue and stiffness matrices due to inertia and + // body forces + // + Array lR(dof,eNoN); + Array3 lK(dof*dof,eNoN,eNoN); + + for (int g = 0; g < lM.nG; g++) { + auto N = lM.N.rcol(g); + double w = lM.w(g) * Jac0 * ht; + + // Acceleration and mass damping at the integration point + // + auto ud = -fb; + + for (int a = 0; a < lM.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 + for (int a = 0; a < lM.eNoN; a++) { + lR(0,a) = lR(0,a) + N(a)*w*ud(0); + lR(1,a) = lR(1,a) + N(a)*w*ud(1); + lR(2,a) = lR(2,a) + N(a)*w*ud(2); + } + + // Local stiffness contribution from mass matrix + for (int b = 0; b < lM.eNoN; b++) { + for (int a = 0; a < lM.eNoN; a++) { + double BtS = w*amd*N(a)*N(b); + lK(0,a,b) = lK(0,a,b) + BtS; + lK(dof+1,a,b) = lK(dof+1,a,b) + BtS; + lK(2*dof+2,a,b) = lK(2*dof+2,a,b) + BtS; + } + } + } + //std::cout << "[shell_cst] " << " " << std::endl; + //std::cout << "[shell_cst] lK: " << lK << std::endl; + //exit(0); + + // Contribution to residue from membrane strain + // + double w = Jac0 * 0.5; + + for (int a = 0; a < lM.eNoN; a++) { + double BtS = Bm(0,0,a)*Sm(0,0) + Bm(1,0,a)*Sm(1,0) + Bm(2,0,a)*Sm(2,0); + lR(0,a) = lR(0,a) + w*BtS; + + BtS = Bm(0,1,a)*Sm(0,0) + Bm(1,1,a)*Sm(1,0) + Bm(2,1,a)*Sm(2,0); + lR(1,a) = lR(1,a) + w*BtS; + + BtS = Bm(0,2,a)*Sm(0,0) + Bm(1,2,a)*Sm(1,0) + Bm(2,2,a)*Sm(2,0); + lR(2,a) = lR(2,a) + w*BtS; + } + + // Contribution to residue from bending strain + // + for (int a = 0; a < eNoN; a++) { + double BtS = 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) + w*BtS; + + BtS = 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) + w*BtS; + + BtS = 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) + w*BtS; + } + + // Contribution to stiffness from membrane-membrane interactions + w = afl * Jac0 * 0.5; + //dmsg << "w: " << w; + + for (int b = 0; b < lM.eNoN; b++) { + for (int a = 0; a < lM.eNoN; a++) { + 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 BtDB = Bm(0,0,a)*D0Bm(0,0,b) + Bm(1,0,a)*D0Bm(1,0,b) + Bm(2,0,a)*D0Bm(2,0,b); + lK(0,a,b) = lK(0,a,b) + w*(BtDB + NxSNx); + + BtDB = Bm(0,0,a)*D0Bm(0,1,b) + Bm(1,0,a)*D0Bm(1,1,b) + Bm(2,0,a)*D0Bm(2,1,b); + lK(1,a,b) = lK(1,a,b) + w*BtDB; + + BtDB = Bm(0,0,a)*D0Bm(0,2,b) + Bm(1,0,a)*D0Bm(1,2,b) + Bm(2,0,a)*D0Bm(2,2,b); + lK(2,a,b) = lK(2,a,b) + w*BtDB; + + BtDB = Bm(0,1,a)*D0Bm(0,0,b) + Bm(1,1,a)*D0Bm(1,0,b) + Bm(2,1,a)*D0Bm(2,0,b); + lK(dof+0,a,b) = lK(dof+0,a,b) + w*BtDB; + + BtDB = Bm(0,1,a)*D0Bm(0,1,b) + Bm(1,1,a)*D0Bm(1,1,b) + Bm(2,1,a)*D0Bm(2,1,b); + lK(dof+1,a,b) = lK(dof+1,a,b) + w*(BtDB + NxSNx); + + BtDB = Bm(0,1,a)*D0Bm(0,2,b) + Bm(1,1,a)*D0Bm(1,2,b) + Bm(2,1,a)*D0Bm(2,2,b); + lK(dof+2,a,b) = lK(dof+2,a,b) + w*BtDB; + + BtDB = Bm(0,2,a)*D0Bm(0,0,b) + Bm(1,2,a)*D0Bm(1,0,b) + Bm(2,2,a)*D0Bm(2,0,b); + lK(2*dof+0,a,b) = lK(2*dof+0,a,b) + w*BtDB; + + BtDB = Bm(0,2,a)*D0Bm(0,1,b) + Bm(1,2,a)*D0Bm(1,1,b) + Bm(2,2,a)*D0Bm(2,1,b); + lK(2*dof+1,a,b) = lK(2*dof+1,a,b) + w*BtDB; + + BtDB = Bm(0,2,a)*D0Bm(0,2,b) + Bm(1,2,a)*D0Bm(1,2,b) + Bm(2,2,a)*D0Bm(2,2,b); + lK(2*dof+2,a,b) = lK(2*dof+2,a,b) + w*(BtDB + NxSNx); + } + } + //std::cout << "[shell_cst] " << " " << std::endl; + //std::cout << "[shell_cst] lK: " << lK << std::endl; + //exit(0); + + // Contribution to stiffness from membrane-bending interactions + // + for (int b = 0; b < eNoN; b++) { + for (int a = 0; a < lM.eNoN; a++) { + double BtDB = Bm(0,0,a)*D1Bb(0,0,b) + Bm(1,0,a)*D1Bb(1,0,b) + Bm(2,0,a)*D1Bb(2,0,b); + lK(0,a,b) = lK(0,a,b) + w*BtDB; + + BtDB = Bm(0,0,a)*D1Bb(0,1,b) + Bm(1,0,a)*D1Bb(1,1,b) + Bm(2,0,a)*D1Bb(2,1,b); + lK(1,a,b) = lK(1,a,b) + w*BtDB; + + BtDB = Bm(0,0,a)*D1Bb(0,2,b) + Bm(1,0,a)*D1Bb(1,2,b) + Bm(2,0,a)*D1Bb(2,2,b); + lK(2,a,b) = lK(2,a,b) + w*BtDB; + + BtDB = Bm(0,1,a)*D1Bb(0,0,b) + Bm(1,1,a)*D1Bb(1,0,b) + Bm(2,1,a)*D1Bb(2,0,b); + lK(dof+0,a,b) = lK(dof+0,a,b) + w*BtDB; + + BtDB = Bm(0,1,a)*D1Bb(0,1,b) + Bm(1,1,a)*D1Bb(1,1,b) + Bm(2,1,a)*D1Bb(2,1,b); + lK(dof+1,a,b) = lK(dof+1,a,b) + w*BtDB; + + BtDB = Bm(0,1,a)*D1Bb(0,2,b) + Bm(1,1,a)*D1Bb(1,2,b) + Bm(2,1,a)*D1Bb(2,2,b); + lK(dof+2,a,b) = lK(dof+2,a,b) + w*BtDB; + + BtDB = Bm(0,2,a)*D1Bb(0,0,b) + Bm(1,2,a)*D1Bb(1,0,b) + Bm(2,2,a)*D1Bb(2,0,b); + lK(2*dof+0,a,b) = lK(2*dof+0,a,b) + w*BtDB; + + BtDB = Bm(0,2,a)*D1Bb(0,1,b) + Bm(1,2,a)*D1Bb(1,1,b) + Bm(2,2,a)*D1Bb(2,1,b); + lK(2*dof+1,a,b) = lK(2*dof+1,a,b) + w*BtDB; + + BtDB = Bm(0,2,a)*D1Bb(0,2,b) + Bm(1,2,a)*D1Bb(1,2,b) + Bm(2,2,a)*D1Bb(2,2,b); + lK(2*dof+2,a,b) = lK(2*dof+2,a,b) + w*BtDB; + } + } + + // Contribution to stiffness from bending-membrane interactions + // + for (int b = 0; b < lM.eNoN; b++) { + for (int a = 0; a < eNoN; a++) { + double BtDB = Bb(0,0,a)*D1Bm(0,0,b) + Bb(1,0,a)*D1Bm(1,0,b) + Bb(2,0,a)*D1Bm(2,0,b); + lK(0,a,b) = lK(0,a,b) + w*BtDB; + + BtDB = Bb(0,0,a)*D1Bm(0,1,b) + Bb(1,0,a)*D1Bm(1,1,b) + Bb(2,0,a)*D1Bm(2,1,b); + lK(1,a,b) = lK(1,a,b) + w*BtDB; + + BtDB = Bb(0,0,a)*D1Bm(0,2,b) + Bb(1,0,a)*D1Bm(1,2,b) + Bb(2,0,a)*D1Bm(2,2,b); + lK(2,a,b) = lK(2,a,b) + w*BtDB; + + BtDB = Bb(0,1,a)*D1Bm(0,0,b) + Bb(1,1,a)*D1Bm(1,0,b) + Bb(2,1,a)*D1Bm(2,0,b); + lK(dof+0,a,b) = lK(dof+0,a,b) + w*BtDB; + + BtDB = Bb(0,1,a)*D1Bm(0,1,b) + Bb(1,1,a)*D1Bm(1,1,b) + Bb(2,1,a)*D1Bm(2,1,b); + lK(dof+1,a,b) = lK(dof+1,a,b) + w*BtDB; + + BtDB = Bb(0,1,a)*D1Bm(0,2,b) + Bb(1,1,a)*D1Bm(1,2,b) + Bb(2,1,a)*D1Bm(2,2,b); + lK(dof+2,a,b) = lK(dof+2,a,b) + w*BtDB; + + BtDB = Bb(0,2,a)*D1Bm(0,0,b) + Bb(1,2,a)*D1Bm(1,0,b) + Bb(2,2,a)*D1Bm(2,0,b); + lK(2*dof+0,a,b) = lK(2*dof+0,a,b) + w*BtDB; + + BtDB = Bb(0,2,a)*D1Bm(0,1,b) + Bb(1,2,a)*D1Bm(1,1,b) + Bb(2,2,a)*D1Bm(2,1,b); + lK(2*dof+1,a,b) = lK(2*dof+1,a,b) + w*BtDB; + + BtDB = Bb(0,2,a)*D1Bm(0,2,b) + Bb(1,2,a)*D1Bm(1,2,b) + Bb(2,2,a)*D1Bm(2,2,b); + lK(2*dof+2,a,b) = lK(2*dof+2,a,b) + w*BtDB; + } + } + + // Contribution to stiffness from bending-bending interactions + // + for (int b = 0; b < eNoN; b++) { + for (int a = 0 ; a < eNoN; a++) { + double BtDB = 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) + w*BtDB; + + BtDB = 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) + w*BtDB; + + BtDB = 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) + w*BtDB; + + BtDB = 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) + w*BtDB; + + BtDB = 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) + w*BtDB; + + BtDB = 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) + w*BtDB; + + BtDB = 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) + w*BtDB; + + BtDB = 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) + w*BtDB; + + BtDB = 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) + w*BtDB; + } + } + + // Global assembly + + //std::cout << "[shell_cst] " << std::endl; + //std::cout << "[shell_cst] lR: " << lR << std::endl; + //std::cout << "[shell_cst] lK: " << lK << std::endl; + //exit(0); + +#ifdef WITH_TRILINOS + if (eq.assmTLS) { + trilinos_doassem_(eNoN, ptr, lK, lR); + } else { +#endif + lhsa_ns::do_assem(com_mod, eNoN, ptr, lK, lR); +#ifdef WITH_TRILINOS + } +#endif +} + //---------- // shell_fp //---------- // void shell_fp(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, - const Array& dl, const Array& xl, const Vector& tfl, Array& lR, Array3& lK) + const Array& dl, const Array& xl, const Array& tfl, Array& lR, Array3& lK) { int nsd = com_mod.nsd; int dof = com_mod.dof; @@ -30,12 +1681,15 @@ void shell_fp(ComMod& com_mod, const int eNoN, const double w, const Vector xc(3,eNoN); double tfn = 0.0; + // [NOTE] This is a hack enabling 'tfl' to be used + // like a vector in the Fortan. + 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(a); + tfn = tfn + N(a)*tfl_data[a]; } double wl = w * tfn; @@ -45,7 +1699,6 @@ void shell_fp(ComMod& com_mod, const int eNoN, const double w, const Vector gCov(3,2); Array gCnv(3,2); nn::gnns(nsd, eNoN, Nx, xc, nV, gCov, gCnv); - //CALL GNNS(eNoN, Nx, xc, nV, gCov, gCnv) // Local residue for (int a = 0; a < eNoN; a++) { @@ -61,24 +1714,126 @@ void shell_fp(ComMod& com_mod, const int eNoN, const double w, const Vector& fNa0, + const double aa_0[2][2], const double aa_x[2][2], const double bb_0[2][2], const double bb_x[2][2], + double& lam3, Array& Sm, Array3& Dm) +{ + using namespace consts; + + #define n_debug_shl_strs_res + #ifdef debug_shl_strs_res + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + #endif + + // Set shell thickness + double ht = lDmn.prop.at(PhysicalProperyType::shell_thickness); + double nu = lDmn.prop.at(PhysicalProperyType::poisson_ratio); + + // Check for incompressibility + bool flag = false; + if (utils::is_zero(nu-0.50)) { + flag = true; + } + //dmsg << "flag: " << flag; + + // Set integration parameters (Gauss coordinates and weights) + double wh[3] = { 5.0/9.0, 5.0/9.0, 8.0/9.0 }; + double xi[3] = { -sqrt(0.60), sqrt(0.60), 0.0 }; + + // Initialize stress-resultants + Sm = 0.0; + Dm = 0.0; + + Vector Sml(3); + Array Dml(3,3); + + // Averaged SQRT(g33) over the thickness + lam3 = 0.0; + + // Gauss integration through shell thickness + // + Array gg_0(2,2), gg_x(2,2); + + for (int g = 0; g < 3; g++) { + //dmsg << "---------- g: " << g+1; + // Local shell thickness coordinate + double xis = 0.50 * ht * xi[g]; + + // Metric coefficients in shell continuum (ref/cur) + for (int i = 0; i < 2; i++) { + for (int j = 0; j < 2; j++) { + gg_0(i,j) = aa_0[i][j] - 2.0*xis*bb_0[i][j]; + gg_x(i,j) = aa_x[i][j] - 2.0*xis*bb_x[i][j]; + } + } + //dmsg << "gg_0: " << gg_0; + //dmsg << "gg_x: " << gg_x; + + // Get 2nd Piola-Kirchhoff and elasticity tensors + // + double g33; + + // For incompressible materials + if (flag) { + mat_models::get_pk2cc_shli(com_mod, lDmn, nFn, fNa0, gg_0, gg_x, g33, Sml, Dml); + + // For compressible materials + } else { + mat_models::get_pk2cc_shlc(com_mod, lDmn, nFn, fNa0, gg_0, gg_x, g33, Sml, Dml); + } + + //dmsg << " " << " "; + //dmsg << "g33: " << g33; + //dmsg << "Sml: " << Sml; + //dmsg << "Dml: " << Dml; + + double wl[3]; + wl[0] = 0.50 * wh[g] * ht; + wl[1] = wl[0] * xis; + wl[2] = wl[1] * xis; + + lam3 = lam3 + wl[0] * sqrt(g33); + + for (int i = 0; i < 3; i++) { + Sm(i,0) = Sm(i,0) + wl[0]*Sml(i); + Sm(i,1) = Sm(i,1) + wl[1]*Sml(i); + + for (int j = 0; j < 3; j++) { + Dm(i,j,0) = Dm(i,j,0) + wl[0]*Dml(i,j); + Dm(i,j,1) = Dm(i,j,1) + wl[1]*Dml(i,j); + Dm(i,j,2) = Dm(i,j,2) + wl[2]*Dml(i,j); + } + } + } + + lam3 = lam3 / ht; + + //dmsg << "Sm: " << Sm; + //dmsg << "Dm: " << Dm; + //dmsg << "lam3: " << lam3; +} }; diff --git a/Code/Source/svFSI/shells.h b/Code/Source/svFSI/shells.h index 511b5bd..8465e08 100644 --- a/Code/Source/svFSI/shells.h +++ b/Code/Source/svFSI/shells.h @@ -7,8 +7,31 @@ namespace shells { +void construct_shell(ComMod& com_mod, const mshType& lM, const Array& Ag, + const Array& Yg, const Array& Dg); + +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); + +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); + +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); + +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& yl, const Array& dl, const Array& xl, + const Array& bfl, const Vector& ptr); + void shell_fp(ComMod& com_mod, const int eNoN, const double w, const Vector& N, const Array& Nx, - const Array& dl, const Array& xl, const Vector& tfl, Array& lR, Array3& lK); + const Array& dl, const Array& xl, const Array& tfl, Array& lR, Array3& lK); + +void shl_strs_res(const ComMod& com_mod, const dmnType& lDmn, const int nFn, const Array& fNa0, + const double aa_0[2][2], const double aa_x[2][2], const double bb_0[2][2], const double bb_x[2][2], + double& lam3, Array& Sm, Array3& Dm); }; diff --git a/Code/Source/svFSI/vtk_xml.cpp b/Code/Source/svFSI/vtk_xml.cpp index cc84757..265d6e5 100644 --- a/Code/Source/svFSI/vtk_xml.cpp +++ b/Code/Source/svFSI/vtk_xml.cpp @@ -881,6 +881,9 @@ void write_vtus(Simulation* simulation, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array& lA, const Array