Add a CMM test case (#117)

* Add a CMM test case. Results are consistent with svFSI output. Prestress examples will be added. * Use a linear solver setting without Trilinos for the cmm inflation test case. * Remove .inp files in the cmm test case. Reduce the number of time steps for the inflation procedure. * Add prestress examples for the cmm test case. --------- Co-authored-by: Martin R. Pfaller <[email protected]>
mrp089 · Oct 9, 2023 · e5d0969 · e5d0969
1 parent aa6a120
commit e5d0969
Show file tree

Hide file tree

Showing 16 changed files with 786 additions and 0 deletions.
diff --git a/tests/cases/cmm_3d_pipe/1-rigid-solution/calcMeanPressTrac.py b/tests/cases/cmm_3d_pipe/1-rigid-solution/calcMeanPressTrac.py
@@ -0,0 +1,151 @@
+# Header
+import numpy as np
+import vtk
+
+#----------------------------------------------------------------------
+def getWallNodes(FileName):
+ modelReader = vtk.vtkXMLPolyDataReader()
+ modelReader.SetFileName(FileName)
+ modelReader.Update()
+
+ model = vtk.vtkPolyData()
+ model = modelReader.GetOutput()
+
+ model_npts = model.GetNumberOfPoints()
+ model_ID = model.GetPointData().GetArray('GlobalNodeID')
+ id_list = np.zeros((model_npts,1), dtype='int32')
+ for i in range(0, model_npts):
+ id_list[i] = int(model_ID.GetTuple1(i))
+ return id_list
+#----------------------------------------------------------------------
+
+#----------------------------------------------------------------------
+def loadVTU(fileName):
+
+ print (" Loading vtu file <--- {}".format(fileName))
+ vtuReader = vtk.vtkXMLUnstructuredGridReader()
+ vtuReader.SetFileName(fileName)
+ vtuReader.Update()
+
+ vtuMesh = vtk.vtkUnstructuredGrid()
+ vtuMesh = vtuReader.GetOutput()
+
+ return vtuMesh
+#----------------------------------------------------------------------
+
+#----------------------------------------------------------------------
+def getSurfaceData(msh, srf_ids, keyWord):
+ msh_npts = msh.GetNumberOfPoints()
+ msh_data = vtk.vtkDoubleArray()
+ msh_data = msh.GetPointData().GetArray(keyWord)
+ num_comp = msh_data.GetNumberOfComponents()
+
+ srf_nno = np.size(srf_ids)
+ if num_comp == 1:
+ srf_data = np.zeros((srf_nno,))
+ else:
+ srf_data = np.zeros((srf_nno,num_comp))
+ for ipt in range(0, srf_nno):
+ h = msh_data.GetTuple(int(srf_ids[ipt])-1)
+ if num_comp == 1:
+ srf_data[ipt] = h[0]
+ else:
+ for j in range(0, num_comp):
+ srf_data[ipt,j] = h[j]
+
+ return srf_data
+#----------------------------------------------------------------------
+
+#----------------------------------------------------------------------
+def writeSrfTraction(h, fout, fwall):
+ modelReader = vtk.vtkXMLPolyDataReader()
+ modelReader.SetFileName(fwall)
+ modelReader.Update()
+
+ model = vtk.vtkPolyData()
+ model = modelReader.GetOutput()
+ model_npts = modelReader.GetNumberOfPoints()
+
+ vtkH = vtk.vtkDoubleArray()
+ vtkH.SetNumberOfComponents(3)
+ vtkH.Allocate(model_npts)
+ vtkH.SetNumberOfTuples(model_npts)
+ vtkH.SetName("Traction")
+ for i in range(0, model_npts):
+ vtkH.SetTuple3(i, h[i,0], h[i,1], h[i,2])
+
+ model.GetPointData().AddArray(vtkH)
+ modelWrite = vtk.vtkXMLPolyDataWriter()
+ modelWrite.SetInputData(model)
+ modelWrite.SetFileName(fout)
+ modelWrite.Write()
+
+ return
+#----------------------------------------------------------------------
+
+#----------------------------------------------------------------------
+def writeSrfPressure(p, fout, fwall):
+ modelReader = vtk.vtkXMLPolyDataReader()
+ modelReader.SetFileName(fwall)
+ modelReader.Update()
+
+ model = vtk.vtkPolyData()
+ model = modelReader.GetOutput()
+ model_npts = modelReader.GetNumberOfPoints()
+
+ vtkP = vtk.vtkDoubleArray()
+ vtkP.SetNumberOfComponents(1)
+ vtkP.Allocate(model_npts)
+ vtkP.SetNumberOfTuples(model_npts)
+ vtkP.SetName("Pressure")
+ for i in range(0, model_npts):
+ vtkP.SetTuple1(i, p[i])
+
+ model.GetPointData().AddArray(vtkP)
+ modelWrite = vtk.vtkXMLPolyDataWriter()
+ modelWrite.SetInputData(model)
+ modelWrite.SetFileName(fout)
+ modelWrite.Write()
+
+ return
+#----------------------------------------------------------------------
+
+if __name__ == '__main__':
+ srcdir = "24-procs"
+ nstart = 600
+ nend = 800
+ nfreq = 100
+ dt = 0.005
+ fwall = "../mesh/mesh-surfaces/lumen_wall.vtp"
+
+ wall_ids = getWallNodes(fwall)
+ wall_nno = np.size(wall_ids)
+ mean_h = np.zeros((wall_nno,3), dtype='float64')
+ mean_P = np.zeros((wall_nno,), dtype='float64')
+
+ nframe = int((nend - nstart)/nfreq) + 1
+ for i in range(0, nframe):
+ ntime = nstart + i*nfreq
+ time = float(ntime)*dt
+ if (ntime < 100):
+ fname = "%s/result_%03d_cpp.vtu" %(srcdir, ntime)
+ else:
+ fname = "%s/result_%d_cpp.vtu" %(srcdir, ntime)
+ print ("Reading file <----- {}".format(fname))
+ vtuMesh = loadVTU(fname)
+ mean_h = mean_h + getSurfaceData(vtuMesh, wall_ids, 'Traction')
+ mean_P = mean_P + getSurfaceData(vtuMesh, wall_ids, 'Pressure')
+
+ mean_h = mean_h / float(nframe)
+ mean_P = mean_P / float(nframe)
+ fout = 'rigid_wall_mean_traction.vtp'
+ print ("Writing traction file ----> {}".format(fout))
+ writeSrfTraction(mean_h, fout, fwall)
+ fout = 'rigid_wall_mean_pressure.vtp'
+ print ("Writing pressure file ----> {}".format(fout))
+ writeSrfPressure(mean_P, fout, fwall)
+
+#=======================================================================
+# EOF
+
+
diff --git a/tests/cases/cmm_3d_pipe/1-rigid-solution/result_800_cpp.vtu b/tests/cases/cmm_3d_pipe/1-rigid-solution/result_800_cpp.vtu
diff --git a/tests/cases/cmm_3d_pipe/1-rigid-solution/rigid_wall_mean_pressure.vtp b/tests/cases/cmm_3d_pipe/1-rigid-solution/rigid_wall_mean_pressure.vtp
diff --git a/tests/cases/cmm_3d_pipe/1-rigid-solution/rigid_wall_mean_traction.vtp b/tests/cases/cmm_3d_pipe/1-rigid-solution/rigid_wall_mean_traction.vtp
diff --git a/tests/cases/cmm_3d_pipe/1-rigid-solution/svFSI.xml b/tests/cases/cmm_3d_pipe/1-rigid-solution/svFSI.xml
@@ -0,0 +1,102 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<svFSIFile version="0.1">
+
+<GeneralSimulationParameters>
+
+ <Continue_previous_simulation> 1 </Continue_previous_simulation>
+ <Number_of_spatial_dimensions> 3 </Number_of_spatial_dimensions> 
+ <Number_of_time_steps> 800 </Number_of_time_steps> 
+ <Time_step_size> 0.005 </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> 1 </Save_results_to_VTK_format> 
+ <Name_prefix_of_saved_VTK_files> result </Name_prefix_of_saved_VTK_files> 
+ <Increment_in_saving_VTK_files> 100 </Increment_in_saving_VTK_files> 
+ <Start_saving_after_time_step> 1 </Start_saving_after_time_step> 
+
+ <Increment_in_saving_restart_files> 100 </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="msh" > 
+
+ <Mesh_file_path> ../../pipe3D_RCR/mesh/mesh-complete.mesh.vtu </Mesh_file_path>
+
+ <Add_face name="lumen_inlet">
+ <Face_file_path> ../../pipe3D_RCR/mesh/mesh-surfaces/lumen_inlet.vtp </Face_file_path>
+ </Add_face>
+
+ <Add_face name="lumen_outlet">
+ <Face_file_path> ../../pipe3D_RCR/mesh/mesh-surfaces/lumen_outlet.vtp </Face_file_path>
+ </Add_face>
+
+ <Add_face name="lumen_wall">
+ <Face_file_path> ../../pipe3D_RCR/mesh/mesh-surfaces/lumen_wall.vtp </Face_file_path>
+ </Add_face> 
+</Add_mesh>
+
+
+<Add_equation type="fluid" > 
+ <Density> 1.06 </Density>
+ <Viscosity model="Constant" >
+ <Value> 0.04 </Value>
+ </Viscosity>
+ <Coupled> true </Coupled>
+ <Min_iterations> 3 </Min_iterations> 
+ <Max_iterations> 10</Max_iterations>
+ <Tolerance> 1e-3 </Tolerance> 
+ <Backflow_stabilization_coefficient> 0.2 </Backflow_stabilization_coefficient> 
+
+ <Output type="Spatial">
+ <Velocity> true </Velocity>
+ <Pressure> true </Pressure>
+ <Traction> true </Traction>
+ <WSS> true </WSS>
+ </Output>
+
+ <LS type="NS" >
+ <Max_iterations> 10 </Max_iterations>
+ <NS_GM_max_iterations> 3 </NS_GM_max_iterations>
+ <NS_CG_max_iterations> 500 </NS_CG_max_iterations>
+ <Tolerance> 1e-3 </Tolerance>
+ <NS_GM_tolerance> 1e-3 </NS_GM_tolerance>
+ <NS_CG_tolerance> 1e-3 </NS_CG_tolerance>
+ <Krylov_space_dimension> 50 </Krylov_space_dimension>
+ </LS>
+
+
+ <Add_BC name="lumen_inlet" > 
+ <Type> Dir </Type>
+ <Time_dependence> Unsteady </Time_dependence>
+ <Temporal_values_file_path> ../lumen_inlet.flow </Temporal_values_file_path>
+ <Profile> Parabolic </Profile>
+ <Impose_flux> true </Impose_flux>
+ </Add_BC> 
+
+ <Add_BC name="lumen_outlet" > 
+ <Type> Neu </Type> 
+ <Time_dependence> RCR </Time_dependence>
+ <RCR_values> 
+ <Proximal_resistance> 121.0 </Proximal_resistance> 
+ <Capacitance> 1.5e-5 </Capacitance> 
+ <Distal_resistance> 1212.0 </Distal_resistance> 
+ <Distal_pressure> 0.0 </Distal_pressure> 
+ <Initial_pressure> 0.0 </Initial_pressure> 
+ </RCR_values>
+ </Add_BC> 
+
+ <Add_BC name="lumen_wall" > 
+ <Type> Dir </Type> 
+ <Time_dependence> Steady </Time_dependence>
+ <Value> 0.0 </Value>
+ </Add_BC> 
+</Add_equation>
+</svFSIFile>
+
+
diff --git a/tests/cases/cmm_3d_pipe/2a-inflate/result_003.vtu b/tests/cases/cmm_3d_pipe/2a-inflate/result_003.vtu
diff --git a/tests/cases/cmm_3d_pipe/2a-inflate/svFSI.xml b/tests/cases/cmm_3d_pipe/2a-inflate/svFSI.xml
@@ -0,0 +1,60 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<svFSIFile version="0.1">
+
+<GeneralSimulationParameters>
+
+ <Continue_previous_simulation> 0 </Continue_previous_simulation>
+ <Number_of_spatial_dimensions> 3 </Number_of_spatial_dimensions> 
+ <Number_of_time_steps> 3 </Number_of_time_steps> 
+ <Time_step_size> 0.001 </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> 1 </Save_results_to_VTK_format> 
+ <Name_prefix_of_saved_VTK_files> result </Name_prefix_of_saved_VTK_files> 
+ <Increment_in_saving_VTK_files> 3 </Increment_in_saving_VTK_files> 
+ <Start_saving_after_time_step> 1 </Start_saving_after_time_step> 
+
+ <Increment_in_saving_restart_files> 3 </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="wall" > 
+ <Set_mesh_as_shell> true </Set_mesh_as_shell>
+ <Mesh_file_path> ../../pipe3D_RCR/mesh/mesh-surfaces/lumen_wall.vtp </Mesh_file_path>
+</Add_mesh>
+
+<Add_equation type="CMM" > 
+ <Coupled> true </Coupled> 
+ <Min_iterations> 3 </Min_iterations> 
+ <Max_iterations> 5 </Max_iterations>
+ <Tolerance> 1e-6 </Tolerance> 
+
+ <Initialize> inflate </Initialize>
+ <Poisson_ratio> 0.5 </Poisson_ratio>
+ <Shell_thickness> 0.2 </Shell_thickness>
+ <Elasticity_modulus> 4000000.0 </Elasticity_modulus>
+
+ <Output type="Spatial">
+ <Displacement> true </Displacement>
+ </Output>
+
+ <LS type="CG" >
+ <Max_iterations> 500 </Max_iterations>
+ <Tolerance> 1e-6 </Tolerance>
+ </LS>
+
+ <Add_BF mesh="wall" > 
+ <Type> traction </Type> 
+ <Time_dependence> spatial </Time_dependence>
+ <Spatial_values_file_path> ../1-rigid-solution/rigid_wall_mean_traction.vtp </Spatial_values_file_path>
+ </Add_BF> 
+</Add_equation>
+</svFSIFile>
+
+
diff --git a/tests/cases/cmm_3d_pipe/2b-prestress/result_003.vtu b/tests/cases/cmm_3d_pipe/2b-prestress/result_003.vtu
diff --git a/tests/cases/cmm_3d_pipe/2b-prestress/svFSI.xml b/tests/cases/cmm_3d_pipe/2b-prestress/svFSI.xml
@@ -0,0 +1,63 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<svFSIFile version="0.1">
+
+<GeneralSimulationParameters>
+
+ <Continue_previous_simulation> 0 </Continue_previous_simulation>
+ <Number_of_spatial_dimensions> 3 </Number_of_spatial_dimensions> 
+ <Number_of_time_steps> 3 </Number_of_time_steps> 
+ <Time_step_size> 0.001 </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> 1 </Save_results_to_VTK_format> 
+ <Name_prefix_of_saved_VTK_files> result </Name_prefix_of_saved_VTK_files> 
+ <Increment_in_saving_VTK_files> 3 </Increment_in_saving_VTK_files> 
+ <Start_saving_after_time_step> 1 </Start_saving_after_time_step> 
+
+ <Increment_in_saving_restart_files> 3 </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="wall" > 
+ <Set_mesh_as_shell> true </Set_mesh_as_shell>
+ <Mesh_file_path> ../../pipe3D_RCR/mesh/mesh-surfaces/lumen_wall.vtp </Mesh_file_path>
+</Add_mesh>
+
+
+<Add_equation type="CMM" > 
+ <Coupled> true </Coupled> 
+ <Min_iterations> 3 </Min_iterations> 
+ <Max_iterations> 5 </Max_iterations>
+ <Tolerance> 1e-6 </Tolerance> 
+ <Prestress> true </Prestress>
+ <Initialize> prestress </Initialize>
+
+ <Poisson_ratio> 0.5 </Poisson_ratio>
+ <Shell_thickness> 0.2 </Shell_thickness>
+ <Elasticity_modulus> 4000000.0 </Elasticity_modulus>
+
+ <Output type="Spatial">
+ <Displacement> true </Displacement>
+ <Stress> true </Stress>
+ </Output>
+
+ <LS type="GMRES" >
+ <Max_iterations> 500 </Max_iterations>
+ <Tolerance> 1e-6 </Tolerance>
+ </LS>
+
+ <Add_BF mesh="wall" > 
+ <Type> traction </Type> 
+ <Time_dependence> spatial </Time_dependence>
+ <Spatial_values_file_path> ../1-rigid-solution/rigid_wall_mean_traction.vtp </Spatial_values_file_path>
+ </Add_BF> 
+</Add_equation>
+</svFSIFile>
+
+