Fsi ustruct

Code/Source/solver/fsi.cpp

@@ -37,6 +37,7 @@
 #include "lhsa.h"
 #include "nn.h"
 #include "sv_struct.h"
+#include "ustruct.h"
 #include "utils.h"
 
 #include <array>
@@ -230,9 +231,9 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar
           break;
 
           case Equation_ustruct:
-            throw std::runtime_error("[construct_fsi] USTRUCT3D_M not implemented");
-            //CALL USTRUCT3D_M(vmsStab, fs(1).eNoN, fs(2).eNoN, nFn, w, Jac, fs(1).N(:,g), fs(2).N(:,g), Nwx, al, yl, 
-            //                 dl, bfl, fN, ya_l, lR, lK, lKd)
+            auto N0 = fs_1[0].N.col(g);
+            auto N1 = fs_1[1].N.col(g);
+            ustruct::ustruct_3d_m(com_mod, cep_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, nFn, w, Jac, N0, N1, Nwx, al, yl, dl, bfl, fN, ya_l, lR, lK, lKd);
           break;
           }
 
@@ -297,8 +298,9 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar
           } break;
 
           case Equation_ustruct:
-            throw std::runtime_error("[construct_fsi] USTRUCT3D_C not implemented");
-            //CALL USTRUCT3D_C(vmsStab, fs(1).eNoN, fs(2).eNoN, w, Jac, fs(1).N(:,g), fs(2).N(:,g), Nwx, Nqx, al, yl, dl, bfl, lR, lK, lKd)
+            auto N0 = fs_2[0].N.col(g);
+            auto N1 = fs_2[1].N.col(g);
+            ustruct::ustruct_3d_c(com_mod, cep_mod, vmsStab, fs_2[0].eNoN, fs_2[1].eNoN, w, Jac, N0, N1, Nwx, Nqx, al, yl, dl, bfl, lR, lK, lKd);
           break;
         }
 
@@ -320,7 +322,15 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar
       }
     } // g: loop
 
-    eq.linear_algebra->assemble(com_mod, eNoN, ptr, lK, lR);
+
+    if (cPhys == Equation_ustruct) 
+    {
+      ustruct::ustruct_do_assem(com_mod, eNoN, ptr, lKd, lK, lR);
+    }
+    else 
+    {
+      eq.linear_algebra->assemble(com_mod, eNoN, ptr, lK, lR);
+    }
 
   } // e: loop
 

Code/Source/solver/main.cpp

@@ -517,7 +517,7 @@ void iterate_solution(Simulation* simulation)
       // Note that this step is done only first iteration. Residual
       // will be 0 for subsequent iterations
       //
-      // Modifies com_mod.Rd.
+      // Modifies com_mod.R.
       //
       #ifdef debug_iterate_solution
       dmsg << "com_mod.sstEq: " << com_mod.sstEq;

tests/cases/fsi/pipe_3d_ustruct/README.md

@@ -0,0 +1,19 @@
+
+# **Problem Description**
+
+Simulate pressure wave propagation in an arterial model using the Arbitrary Lagrangian-Eulerian method [1]. The problem set-up is as follows.
+
+<p align="center">
+   <img src="./configuration.png" width="600">
+</p>
+
+And the results are
+
+<p align="center">
+   <img src="./results.gif" width="600">
+</p>
+
+
+## References
+
+1.  Liu, Ju, and Alison L. Marsden.  A Unified Continuum and Variational Multiscale Formulation for Fluids, Solids, and Fluid Structure Interaction.  *Computer Methods in Applied Mechanics and Engineering* 337 (August 2018): 549 97. https://doi.org/10.1016/j.cma.2018.03.045.

tests/cases/fsi/pipe_3d_ustruct/solver.xml

@@ -0,0 +1,156 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<svMultiPhysicsFile version="0.1">
+
+<GeneralSimulationParameters>
+  <Continue_previous_simulation> 0 </Continue_previous_simulation>
+  <Number_of_spatial_dimensions> 3 </Number_of_spatial_dimensions> 
+  <Number_of_time_steps> 5 </Number_of_time_steps> 
+  <Time_step_size> 1e-4 </Time_step_size> 
+  <Spectral_radius_of_infinite_time_step> 0.50 </Spectral_radius_of_infinite_time_step> 
+  <Searched_file_name_to_trigger_stop> STOP_SIM </Searched_file_name_to_trigger_stop> 
+  <Save_results_to_VTK_format> true </Save_results_to_VTK_format> 
+  <Name_prefix_of_saved_VTK_files> result </Name_prefix_of_saved_VTK_files> 
+  <Increment_in_saving_VTK_files> 5 </Increment_in_saving_VTK_files> 
+  <Start_saving_after_time_step> 1 </Start_saving_after_time_step> 
+  <Increment_in_saving_restart_files> 1 </Increment_in_saving_restart_files> 
+  <Convert_BIN_to_VTK_format> 0 </Convert_BIN_to_VTK_format> 
+  <Verbose> 1 </Verbose> 
+  <Warning> 0 </Warning> 
+  <Debug> 0 </Debug> 
+</GeneralSimulationParameters>
+
+<Add_mesh name="lumen" > 
+  <Mesh_file_path> mesh/fluid/mesh-complete.mesh.vtu  </Mesh_file_path>
+  <Add_face name="lumen_inlet">
+      <Face_file_path> mesh/fluid/mesh-surfaces/start.vtp </Face_file_path>
+  </Add_face>
+  <Add_face name="lumen_outlet">
+      <Face_file_path> mesh/fluid/mesh-surfaces/end.vtp </Face_file_path>
+  </Add_face>
+  <Add_face name="lumen_wall">
+      <Face_file_path> mesh/fluid/mesh-surfaces/interface.vtp </Face_file_path>
+  </Add_face>
+  <Domain> 0 </Domain>
+
+</Add_mesh>
+
+<Add_mesh name="wall" >
+  <Mesh_file_path> mesh/solid/mesh-complete.mesh.vtu  </Mesh_file_path>
+  <Add_face name="wall_inlet">
+      <Face_file_path> mesh/solid/mesh-surfaces/start.vtp </Face_file_path>
+  </Add_face>
+  <Add_face name="wall_outlet">
+      <Face_file_path> mesh/solid/mesh-surfaces/end.vtp </Face_file_path>
+  </Add_face>
+  <Add_face name="wall_inner">
+      <Face_file_path> mesh/solid/mesh-surfaces/interface.vtp </Face_file_path>
+  </Add_face>
+  <Add_face name="wall_outer">
+      <Face_file_path> mesh/solid/mesh-surfaces/outside.vtp </Face_file_path>
+  </Add_face>
+  <Domain> 1 </Domain>
+</Add_mesh>
+
+<Add_projection name="wall_inner" >
+   <Project_from_face> lumen_wall </Project_from_face>
+</Add_projection> 
+
+<Add_equation type="FSI" > 
+   <Coupled> true </Coupled>
+   <Min_iterations> 1 </Min_iterations>  
+   <Max_iterations> 7 </Max_iterations> 
+   <Tolerance> 1e-12 </Tolerance> 
+
+   <Domain id="0" >
+      <Equation> fluid </Equation> 
+      <Density> 1.0 </Density> 
+      <Viscosity model="Constant" >
+         <Value> 0.04 </Value>
+      </Viscosity>
+      <Backflow_stabilization_coefficient> 0.2 </Backflow_stabilization_coefficient> 
+   </Domain>
+
+   <Domain id="1" >
+      <Equation> ustruct </Equation> 
+      <Constitutive_model type="neoHookean"> </Constitutive_model> 
+      <Dilational_penalty_model> M94 </Dilational_penalty_model> 
+      <Density> 1.0 </Density> 
+      <Elasticity_modulus> 1.0e7 </Elasticity_modulus> 
+      <Poisson_ratio> 0.3 </Poisson_ratio> 
+   </Domain>
+
+   <LS type="GMRES" >
+      <Linear_algebra type="fsils" >
+         <Preconditioner> fsils </Preconditioner>
+      </Linear_algebra>
+      <Tolerance> 1e-12 </Tolerance>
+      <Max_iterations> 100 </Max_iterations> 
+      <Krylov_space_dimension> 50 </Krylov_space_dimension>
+   </LS>
+
+   <Output type="Spatial" >
+     <Displacement> true </Displacement>
+     <Velocity> true </Velocity>
+     <Pressure> true </Pressure>
+     <VonMises_stress> true </VonMises_stress>
+   </Output>
+
+   <Output type="Alias" >
+       <Displacement> FS_Displacement </Displacement>
+   </Output>
+
+   <Add_BC name="lumen_inlet" > 
+      <Type> Neu </Type> 
+      <Value> 5.0e4 </Value> 
+   </Add_BC> 
+
+   <Add_BC name="wall_inlet" > 
+      <Type> Dir </Type> 
+      <Value> 0.0 </Value> 
+      <Impose_on_state_variable_integral> true </Impose_on_state_variable_integral> 
+      <Zero_out_perimeter> false </Zero_out_perimeter> 
+      <Effective_direction> (0, 0, 1) </Effective_direction> 
+   </Add_BC> 
+
+   <Add_BC name="wall_outlet" >
+      <Type> Dir </Type> 
+      <Value> 0.0 </Value> 
+      <Impose_on_state_variable_integral> true </Impose_on_state_variable_integral>
+      <Zero_out_perimeter> false </Zero_out_perimeter> 
+      <Effective_direction> (0, 0, 1 ) </Effective_direction> 
+   </Add_BC> 
+
+</Add_equation>
+
+
+<Add_equation type="mesh" >
+   <Coupled> true </Coupled>
+   <Min_iterations> 1 </Min_iterations>
+   <Max_iterations> 7 </Max_iterations>
+   <Tolerance> 1e-12 </Tolerance>
+   <Poisson_ratio> 0.3 </Poisson_ratio> 
+
+   <LS type="CG" >
+      <Linear_algebra type="fsils" >
+         <Preconditioner> fsils </Preconditioner>
+      </Linear_algebra>
+      <Tolerance> 1e-12 </Tolerance>
+   </LS>
+
+   <Output type="Spatial" >
+     <Displacement> true </Displacement>
+   </Output>
+
+   <Add_BC name="lumen_inlet" > 
+      <Type> Dir </Type> 
+      <Value> 0.0 </Value> 
+   </Add_BC> 
+
+   <Add_BC name="lumen_outlet" > 
+      <Type> Dir </Type> 
+      <Value> 0.0 </Value> 
+   </Add_BC> 
+
+</Add_equation>
+
+</svMultiPhysicsFile>

tests/cases/fsi/pipe_RCR_3d_ustruct/README.md

@@ -0,0 +1,3 @@
+
+# **Problem Description**
+Same as pipe_3d, but uses RCR boundary condition on lumen_outlet.