fsi_tissue: 在已有的fsi程序中引入了用于生成局部纤维的局部坐标系功能。

examples/fsi/CMakeLists.txt

@@ -62,6 +62,7 @@ SET( perigee_preprocess_lib_src
   ${perigee_source}/Mesh/Global_Part_Reload.cpp
   ${perigee_source}/Mesh/Part_FEM.cpp
   ${perigee_source}/Mesh/Part_FEM_FSI.cpp
+  ${perigee_source}/Mesh/Part_FEM_FSI_Tissue.cpp
   ${perigee_source}/Element/FE_Tools.cpp
   ${perigee_source}/Element/FEAElement_Triangle3_3D_der0.cpp
   ${perigee_source}/Element/FEAElement_Triangle6_3D_der0.cpp
@@ -103,7 +104,10 @@ SET( perigee_analysis_lib_src
   ${perigee_source}/Element/FEAElement_Hex8.cpp
   ${perigee_source}/Model/MaterialModel_NeoHookean_M94_Mixed.cpp
   ${perigee_source}/Model/MaterialModel_NeoHookean_Incompressible_Mixed.cpp
+  ${perigee_source}/Model/MaterialModel_GOH06_Incompressible_Mixed.cpp
+  ${perigee_source}/Model/MaterialModel_GOH06_ST91_Mixed.cpp
   ${perigee_source}/Model/Tissue_prestress.cpp
+  ${perigee_source}/Model/Tissue_property.cpp
   ${perigee_source}/Model/GenBC_RCR.cpp
   ${perigee_source}/Model/GenBC_Resistance.cpp
   ${perigee_source}/Model/GenBC_Inductance.cpp

examples/fsi/fsi_driver.cpp

@@ -27,6 +27,8 @@
 #include "GenBC_Pressure.hpp"
 #include "MaterialModel_NeoHookean_M94_Mixed.hpp"
 #include "MaterialModel_NeoHookean_Incompressible_Mixed.hpp"
+#include "MaterialModel_GOH06_ST91_Mixed.hpp"
+#include "MaterialModel_GOH06_Incompressible_Mixed.hpp"
 #include "PLocAssem_2x2Block_ALE_VMS_NS_GenAlpha.hpp"
 #include "PLocAssem_2x2Block_VMS_Incompressible.hpp"
 #include "PLocAssem_2x2Block_VMS_Hyperelasticity.hpp"
@@ -35,7 +37,6 @@
 #include "PGAssem_FSI.hpp"
 #include "PGAssem_Mesh.hpp"
 #include "PTime_FSI_Solver.hpp"
-
 #include "PDNTimeStep.hpp"
 #include "PETSc_Tools.hpp"
 #include "APart_Node_FSI.hpp"
@@ -59,6 +60,14 @@ int main(int argc, char *argv[])
   double solid_density = 1.0;
   double solid_E = 2.0e6;
   double solid_nu = 0.5;
+  double solid_mu = 6.67e5;
+  double solid_f1the = 40.02;
+  double solid_f1phi = 0.0;
+  double solid_f2the = -40.02;
+  double solid_f2phi = 0.0;
+  double solid_fk1 = 9.966e5;
+  double solid_fk2 = 524.6;
+  double solid_fkd = 0.333;
 
   // mesh motion elasticity solver parameters
   double mesh_E  = 1.0;
@@ -159,6 +168,14 @@ int main(int argc, char *argv[])
   SYS_T::GetOptionReal(  "-sl_density",        solid_density);
   SYS_T::GetOptionReal(  "-sl_E",              solid_E);
   SYS_T::GetOptionReal(  "-sl_nu",             solid_nu);
+  SYS_T::GetOptionReal(  "-sl_mu",             solid_mu);
+  SYS_T::GetOptionReal(  "-sl_f1the",          solid_f1the);
+  SYS_T::GetOptionReal(  "-sl_f1phi",          solid_f1phi);
+  SYS_T::GetOptionReal(  "-sl_f2the",          solid_f2the);
+  SYS_T::GetOptionReal(  "-sl_f2phi",          solid_f2phi);
+  SYS_T::GetOptionReal(  "-sl_fk1",            solid_fk1);
+  SYS_T::GetOptionReal(  "-sl_fk2",            solid_fk2);
+  SYS_T::GetOptionReal(  "-sl_fkd",            solid_fkd);
   SYS_T::GetOptionReal(  "-mesh_E",            mesh_E);
   SYS_T::GetOptionReal(  "-mesh_nu",           mesh_nu);
   SYS_T::GetOptionInt(   "-inflow_type",       inflow_type);
@@ -198,6 +215,14 @@ int main(int argc, char *argv[])
   SYS_T::cmdPrint("-sl_density:", solid_density);
   SYS_T::cmdPrint("-sl_E:", solid_E);
   SYS_T::cmdPrint("-sl_nu:", solid_nu);
+  SYS_T::cmdPrint("-sl_mu:", solid_mu);
+  SYS_T::cmdPrint("-sl_f1the:", solid_f1the);
+  SYS_T::cmdPrint("-sl_f1phi:", solid_f1phi);
+  SYS_T::cmdPrint("-sl_f2the:", solid_f2the);
+  SYS_T::cmdPrint("-sl_f2phi:", solid_f2phi);
+  SYS_T::cmdPrint("-sl_fk1:", solid_fk1);
+  SYS_T::cmdPrint("-sl_fk2:", solid_fk2);
+  SYS_T::cmdPrint("-sl_fkd:", solid_fkd);
   SYS_T::cmdPrint("-mesh_E:", mesh_E);
   SYS_T::cmdPrint("-mesh_nu:", mesh_nu);
 
@@ -331,6 +356,8 @@ int main(int argc, char *argv[])
 
   Tissue_prestress * ps_data = new Tissue_prestress(locElem, nqp_vol, rank, is_load_ps, "./ps_data/prestress");
 
+  Tissue_property * tp_data = new Tissue_property(part_v_file, rank);
+
   // Group APart_Node and ALocal_NBC into a vector
   std::vector<APart_Node *> pNode_list { pNode_v, pNode_p };
 
@@ -451,14 +478,20 @@ int main(int argc, char *argv[])
 
   if( solid_nu == 0.5 )
   {
-    matmodel = new MaterialModel_NeoHookean_Incompressible_Mixed( solid_density, solid_E );
+    matmodel = new MaterialModel_GOH06_Incompressible_Mixed( solid_density, solid_mu,
+        solid_f1the, solid_f1phi, solid_f2the, solid_f2phi, solid_fk1, solid_fk2, solid_fkd );
+
+    //matmodel = new MaterialModel_NeoHookean_Incompressible_Mixed( solid_density, solid_E );
 
     locAssem_solid_ptr = new PLocAssem_2x2Block_VMS_Incompressible(
         matmodel, tm_galpha_ptr, elementv -> get_nLocBas(), elements->get_nLocBas() );
   }
   else
   {
-    matmodel = new MaterialModel_NeoHookean_M94_Mixed( solid_density, solid_E, solid_nu );
+    matmodel = new MaterialModel_GOH06_ST91_Mixed( solid_density, solid_E, solid_nu,
+        solid_f1the, solid_f1phi, solid_f2the, solid_f2phi, solid_fk1, solid_fk2, solid_fkd );
+
+    //matmodel = new MaterialModel_NeoHookean_M94_Mixed( solid_density, solid_E, solid_nu );
 
     locAssem_solid_ptr = new PLocAssem_2x2Block_VMS_Hyperelasticity(
         matmodel, tm_galpha_ptr, elementv -> get_nLocBas(), elements->get_nLocBas() );
@@ -593,7 +626,7 @@ int main(int argc, char *argv[])
 
     gloAssem_ptr->Assem_mass_residual( disp, velo, pres, locElem, locAssem_fluid_ptr,
         locAssem_solid_ptr, elementv, elements, quadv, quads, locIEN_v, locIEN_p, 
-        fNode, locnbc_v, locnbc_p, locebc_v, ps_data );
+        fNode, locnbc_v, locnbc_p, locebc_v, ps_data, tp_data );
 
     Vec proj_vp, proj_v, proj_p;
     VecDuplicate( gloAssem_ptr->G, &proj_vp );
@@ -746,7 +779,7 @@ int main(int argc, char *argv[])
       tm_galpha_ptr, timeinfo, inflow_rate_ptr, locElem, locIEN_v, locIEN_p, 
       pNode_v, pNode_p, fNode, locnbc_v, locnbc_p, locinfnbc, mesh_locnbc, 
       locebc_v, locebc_p, mesh_locebc, 
-      gbc, pmat, mmat, elementv, elements, quadv, quads, ps_data,
+      gbc, pmat, mmat, elementv, elements, quadv, quads, ps_data, tp_data,
       locAssem_fluid_ptr, locAssem_solid_ptr, locAssem_mesh_ptr,
       gloAssem_ptr, gloAssem_mesh_ptr, lsolver, mesh_lsolver, nsolver);
 
@@ -771,7 +804,7 @@ int main(int argc, char *argv[])
   delete GMIptr; delete PartBasic; delete locElem; delete fNode; delete pNode_v; delete pNode_p;
   delete locinfnbc; delete locnbc_v; delete locnbc_p; delete mesh_locnbc; 
   delete locebc_v; delete locebc_p; delete mesh_locebc; 
-  delete locIEN_v; delete locIEN_p; delete ps_data;
+  delete locIEN_v; delete locIEN_p; delete ps_data; delete tp_data;
   PetscFinalize();
   return EXIT_SUCCESS;
 }

examples/fsi/include/PGAssem_FSI.hpp

@@ -65,7 +65,8 @@ class PGAssem_FSI : public IPGAssem
         const ALocal_NBC * const &nbc_v,
         const ALocal_NBC * const &nbc_p,
         const ALocal_EBC * const &ebc_part,
-        const Tissue_prestress * const &ps_ptr );
+        const Tissue_prestress * const &ps_ptr,
+        const Tissue_property * const &tp_ptr );
 
     virtual void Assem_Residual(
         const double &curr_time, const double &dt,
@@ -92,7 +93,8 @@ class PGAssem_FSI : public IPGAssem
         const ALocal_NBC * const &nbc_p,
         const ALocal_EBC * const &ebc_part,
         const IGenBC * const &gbc,
-        const Tissue_prestress * const &ps_ptr );
+        const Tissue_prestress * const &ps_ptr,
+        const Tissue_property * const &tp_ptr );
 
     virtual void Assem_Tangent_Residual(
         const double &curr_time, const double &dt,
@@ -119,7 +121,8 @@ class PGAssem_FSI : public IPGAssem
         const ALocal_NBC * const &nbc_p,
         const ALocal_EBC * const &ebc_part,
         const IGenBC * const &gbc,
-        const Tissue_prestress * const &ps_ptr );
+        const Tissue_prestress * const &ps_ptr,
+	const Tissue_property * const &tp_ptr );
 
     // Assembly routine for the surface integrals for flow rates
     // and averaged pressure

examples/fsi/include/PGAssem_Wall_Prestress.hpp

@@ -57,7 +57,8 @@ class PGAssem_Wall_Prestress : public IPGAssem
         const ALocal_NBC * const &nbc_p,
         const ALocal_EBC * const &ebc_v,
         const ALocal_EBC * const &ebc_p,
-        const Tissue_prestress * const &ps_ptr );
+        const Tissue_prestress * const &ps_ptr,
+        const Tissue_property * const &tp_ptr );
 
     virtual void Assem_Tangent_Residual(
         const double &curr_time,
@@ -81,7 +82,8 @@ class PGAssem_Wall_Prestress : public IPGAssem
         const ALocal_NBC * const &nbc_p,
         const ALocal_EBC * const &ebc_v,
         const ALocal_EBC * const &ebc_p,
-        const Tissue_prestress * const &ps_ptr );
+        const Tissue_prestress * const &ps_ptr,
+        const Tissue_property * const &tp_ptr );
 
     virtual void Update_Wall_Prestress(
         const PDNSolution * const &disp,
@@ -93,7 +95,8 @@ class PGAssem_Wall_Prestress : public IPGAssem
         const ALocal_IEN * const &lien_v,
         const ALocal_IEN * const &lien_p,
         const FEANode * const &fnode_ptr,
-        Tissue_prestress * const &ps_ptr ) const;
+        Tissue_prestress * const &ps_ptr,
+	const Tissue_property * const &tp_ptr ) const;
 
   private:
     const int nLocBas, snLocBas, num_ebc, nlgn_v, nlgn_p;

examples/fsi/include/PLocAssem_2x2Block_VMS_Hyperelasticity.hpp

@@ -59,7 +59,10 @@ class PLocAssem_2x2Block_VMS_Hyperelasticity : public IPLocAssem_2x2Block
         const double * const &eleCtrlPts_y,
         const double * const &eleCtrlPts_z,
         const double * const &qua_prestress,
-        const IQuadPts * const &quad );
+        const IQuadPts * const &quad,
+        const std::vector<Vector_3> &eleBasis_r,
+        const std::vector<Vector_3> &eleBasis_c,
+        const std::vector<Vector_3> &eleBasis_l );
 
     virtual void Assem_Tangent_Residual(
         const double &time, const double &dt,
@@ -74,7 +77,10 @@ class PLocAssem_2x2Block_VMS_Hyperelasticity : public IPLocAssem_2x2Block
         const double * const &eleCtrlPts_y,
         const double * const &eleCtrlPts_z,
         const double * const &qua_prestress,
-        const IQuadPts * const &quad );
+        const IQuadPts * const &quad,
+        const std::vector<Vector_3> &eleBasis_r,
+        const std::vector<Vector_3> &eleBasis_c,
+        const std::vector<Vector_3> &eleBasis_l );
 
     virtual void Assem_Mass_Residual(
         const double * const &disp,
@@ -85,7 +91,10 @@ class PLocAssem_2x2Block_VMS_Hyperelasticity : public IPLocAssem_2x2Block
         const double * const &eleCtrlPts_y,
         const double * const &eleCtrlPts_z,
         const double * const &qua_prestress,
-        const IQuadPts * const &quad );
+        const IQuadPts * const &quad, 
+        const std::vector<Vector_3> &eleBasis_r,
+        const std::vector<Vector_3> &eleBasis_c,
+        const std::vector<Vector_3> &eleBasis_l );
 
     virtual void Assem_Residual_EBC(
         const int &ebc_id,
@@ -116,15 +125,18 @@ class PLocAssem_2x2Block_VMS_Hyperelasticity : public IPLocAssem_2x2Block
         const double * const &eleCtrlPts_x,
         const double * const &eleCtrlPts_y,
         const double * const &eleCtrlPts_z,
-        const IQuadPts * const &quad ) const;
+        const IQuadPts * const &quad,
+        const std::vector<Vector_3> &eleBasis_r,
+        const std::vector<Vector_3> &eleBasis_c,
+        const std::vector<Vector_3> &eleBasis_l ) const;
 
   private:
     const double rho0, alpha_f, alpha_m, gamma;
 
     const int nLocBas, snLocBas, vec_size_0, vec_size_1, sur_size_0;
 
     // useful tensors for the material model
-    const IMaterialModel * const matmodel;
+    IMaterialModel * const matmodel;
 
     void print_info() const;
 

examples/fsi/include/PLocAssem_2x2Block_VMS_Incompressible.hpp

@@ -59,7 +59,10 @@ class PLocAssem_2x2Block_VMS_Incompressible : public IPLocAssem_2x2Block
         const double * const &eleCtrlPts_y,
         const double * const &eleCtrlPts_z,
         const double * const &qua_prestress,
-        const IQuadPts * const &quad );
+        const IQuadPts * const &quad,
+        const std::vector<Vector_3> &eleBasis_r,
+        const std::vector<Vector_3> &eleBasis_c,
+        const std::vector<Vector_3> &eleBasis_l );
 
     virtual void Assem_Tangent_Residual(
         const double &time, const double &dt,
@@ -74,7 +77,10 @@ class PLocAssem_2x2Block_VMS_Incompressible : public IPLocAssem_2x2Block
         const double * const &eleCtrlPts_y,
         const double * const &eleCtrlPts_z,
         const double * const &qua_prestress,
-        const IQuadPts * const &quad );
+        const IQuadPts * const &quad,
+        const std::vector<Vector_3> &eleBasis_r,
+        const std::vector<Vector_3> &eleBasis_c,
+        const std::vector<Vector_3> &eleBasis_l );
 
     virtual void Assem_Mass_Residual(
         const double * const &disp,
@@ -85,7 +91,10 @@ class PLocAssem_2x2Block_VMS_Incompressible : public IPLocAssem_2x2Block
         const double * const &eleCtrlPts_y,
         const double * const &eleCtrlPts_z,
         const double * const &qua_prestress,
-        const IQuadPts * const &quad );
+        const IQuadPts * const &quad,
+       	const std::vector<Vector_3> &eleBasis_r,
+        const std::vector<Vector_3> &eleBasis_c,
+        const std::vector<Vector_3> &eleBasis_l );
 
     virtual void Assem_Residual_EBC(
         const int &ebc_id,
@@ -116,7 +125,10 @@ class PLocAssem_2x2Block_VMS_Incompressible : public IPLocAssem_2x2Block
         const double * const &eleCtrlPts_x,
         const double * const &eleCtrlPts_y,
         const double * const &eleCtrlPts_z,
-        const IQuadPts * const &quad ) const;
+        const IQuadPts * const &quad,
+        const std::vector<Vector_3> &eleBasis_r,
+        const std::vector<Vector_3> &eleBasis_c,
+        const std::vector<Vector_3> &eleBasis_l ) const;
 
   private:
     const double rho0, alpha_f, alpha_m, gamma;
@@ -126,7 +138,7 @@ class PLocAssem_2x2Block_VMS_Incompressible : public IPLocAssem_2x2Block
     const int nLocBas, snLocBas, vec_size_0, vec_size_1, sur_size_0;
 
     // useful tensors for the material model
-    const IMaterialModel * const matmodel;
+    IMaterialModel * const matmodel;
 
     void print_info() const;
 

examples/fsi/include/PNonlinear_FSI_Solver.hpp

@@ -64,6 +64,7 @@ class PNonlinear_FSI_Solver
         const IQuadPts * const &quad_v,
         const IQuadPts * const &quad_s,
         const Tissue_prestress * const &ps_ptr,
+	const Tissue_property * const &tp_ptr,
         IPLocAssem_2x2Block * const &lassem_fluid_ptr,
         IPLocAssem_2x2Block * const &lassem_solid_ptr,
         IPLocAssem * const &lassem_mesh_ptr,
@@ -109,6 +110,7 @@ class PNonlinear_FSI_Solver
         const IQuadPts * const &quad_v,
         const IQuadPts * const &quad_s,
         Tissue_prestress * const &ps_ptr,
+	const Tissue_property * const &tp_ptr,
         IPLocAssem_2x2Block * const &lassem_solid_ptr,
         IPGAssem * const &gassem_ptr,
         PLinear_Solver_PETSc * const &lsolver_ptr,

examples/fsi/include/PTime_FSI_Solver.hpp

@@ -55,7 +55,8 @@ class PTime_FSI_Solver
         FEAElement * const &elements,
         const IQuadPts * const &quad_v,
         const IQuadPts * const &quad_s,
-        const Tissue_prestress * const &ps_ptr,
+        const Tissue_prestress * const &ps_ptr, 
+	const Tissue_property * const &tp_ptr,
         IPLocAssem_2x2Block * const &lassem_fluid_ptr,
         IPLocAssem_2x2Block * const &lassem_solid_ptr,
         IPLocAssem * const &lassem_mesh_ptr,
@@ -94,6 +95,7 @@ class PTime_FSI_Solver
         const IQuadPts * const &quad_v,
         const IQuadPts * const &quad_s,
         Tissue_prestress * const &ps_ptr,
+	const Tissue_property * const &tp_ptr,
         IPLocAssem_2x2Block * const &lassem_solid_ptr,
         IPGAssem * const &gassem_ptr,
         PLinear_Solver_PETSc * const &lsolver_ptr,

examples/fsi/preprocess_fsi.cpp

@@ -14,7 +14,7 @@
 #include "Global_Part_METIS.hpp"
 #include "Global_Part_Serial.hpp"
 #include "Global_Part_Reload.hpp"
-#include "Part_FEM_FSI.hpp"
+#include "Part_FEM_FSI_Tissue.hpp"
 #include "NodalBC.hpp"
 #include "NodalBC_3D_FSI.hpp"
 #include "NodalBC_3D_inflow.hpp"
@@ -394,6 +394,22 @@ int main( int argc, char * argv[] )
   std::cout<<"Solid domain: "<<v_node_s.size()<<" nodes.\n";
   std::cout<<"Fluid-Solid interface: "<<nFunc_interface<<" nodes.\n";
 
+  // --------------------------------------------------------------------------
+  // Read the geometry file for the solid domain, generate the list of direction
+  // basis vectors of the nodes. The list includes radial, longitudinal, and
+  // circumferential basis, denoting by r, l, and c, respectively.
+  const std::vector<int> solid_node_id = VTK_T::read_int_PointData(geo_s_file, "GlobalNodeID");
+  const std::vector<Vector_3> basis_r  = VTK_T::read_Vector_3_PointData(geo_s_file, "radial_basis");
+  const std::vector<Vector_3> basis_c  = VTK_T::read_Vector_3_PointData(geo_s_file, "circumferential_basis");
+  const std::vector<Vector_3> basis_l  = VTK_T::read_Vector_3_PointData(geo_s_file, "longitudinal_basis");
+
+  SYS_T::print_fatal_if(v_node_s != solid_node_id, "ERROR: GlobalNodeID for solid geometry file is not equal to the whole FSI domain.");
+  SYS_T::print_fatal_if(solid_node_id.size() != basis_r.size(), "ERROR: radial_basis is not matched.");
+  SYS_T::print_fatal_if(solid_node_id.size() != basis_c.size(), "ERROR: circumferential_basis is not matched.");
+  SYS_T::print_fatal_if(solid_node_id.size() != basis_l.size(), "ERROR: longitudinal_basis is not matched.");
+  std::cout<<"=== Direction basis vectors generated.\n";
+  // --------------------------------------------------------------------------
+
   std::vector<IIEN const *> ienlist;
   ienlist.push_back(IEN_p); ienlist.push_back(IEN_v);
 
@@ -566,9 +582,9 @@ int main( int argc, char * argv[] )
 
     part_p -> write( part_file_p );
 
-    IPart * part_v = new Part_FEM_FSI( mesh_v, global_part, mnindex_v, IEN_v,
-        ctrlPts, phy_tag, v_node_f, v_node_s, proc_rank, cpu_size, elemType, 
-        start_idx_v[proc_rank], { 1, dof_fields[1], true, "velocity"} );
+    IPart * part_v = new Part_FEM_FSI_Tissue( mesh_v, global_part, mnindex_v, IEN_v,
+        ctrlPts, phy_tag, v_node_f, v_node_s, basis_r, basis_c, basis_l,
+	proc_rank, cpu_size, elemType, start_idx_v[proc_rank], { 1, dof_fields[1], true, "velocity"} );
 
     part_v -> print_part_loadbalance_edgecut();