可视化软件 Paraview 处理 Peridynamics(近场动力学)计算结果
[1] 什么是 peridynamics ? 近场动力学，一种新起的非局部化模型；可用无网格法和有限元法求解，得到的结果需要后处理
[2] Paraview，一款可视化软件，可以将无网格法和有限元法求解得到的结果可视化，得到三维结构。
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
本文关注：用无网格法求解peridynamics，采用fortran语言在VS2013上编写，输出paraview所需要的.vtk文件，最后用paraview可视化
为了达到上述目的，需要完成以下几点：
1. 了解Paraview所需要的 .vtk文件是什么样的格式
2. 用fortran编写peridynamics模型的求解程序，将计算得到的坐标信息，位移，速度，损伤等以vtk格式保存
3. 用paraview 看看生成的vtk格式文件
文末附生成VTK文件的子程序，供大家参考使用。
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
本文算例：冲击有预先裂缝的薄板，来自试验
Kalthoff JF, Winkler S (1988) Failure mode transition at high rates of shear loading. In: Chiem
CY, Kunze H-D, Meyer LW (eds) Impact loading and dynamic behavior of materials, vol
1. DGM Informationsgesellschaft Verlag, Oberursel, pp 185–195
这是文献里面的结果：
本文模拟得到的损伤云图：
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
1. 了解Paraview所需要的 .vtk文件是什么样的格式
参考VTK格式官网：https://lorensen.github.io/VTKExamples/site/VTKFileFormats/
本文只关注无网格，也就是只有粒子，没有cell等信息的一种，总结如下：
# vtk DataFile Version 3.0 ----------- 第一行，表明vtk文件的版本，可以是3.0，4.0等；
title_Impact ----------- 第二行，标题，随便取，自己好记
ASCII ----------- 第三行，数字格式，可以是ASCII或者二进制，我喜欢ASCII直白呀；二进制还得去转换，麻烦
DATASET UNSTRUCTURED_grid ------------第四行，表明几何拓扑，可以去官网对应看一下，有总共五种；这种对于无网格的很好用
POINTS n float -------------第五行，点，点的数量，点坐标的格式，float浮点数
1x，1y, 1z
ix, iy, iz ------------- 这里是点的坐标，从第一个到第n个
nx，ny，nz
Point_data n --------------- 这一行一定不能缺失，而且n要跟上头n一样，都等于节点的数量
SCALARS DAMAGE FLOAT 1 -------------- 记录点上的标量信息，比如损伤；SCALARS DAMAGE FLOAT 1 分别表示标量，标量的名字，自己取；
标量的数字格式，float浮点数，后面的1不能少，不知道干啥的；
LOOKUP_TABLE DEFAULT ------------- 如果记录的是标量，这行一定要有，用来将标量映射为颜色云图；用默认的就行；
d1
di --------------- n个损伤标量信息，相应于n个节点
dn
VECTORS DISPLACEMENT FLOAT ----------- 记录点上的矢量信息，这里是位移；
disp1x，disp1y, disp1z
dispix, dispiy, dispiz ----------记录的n对位移矢量
dispnx, dispny, dispnz
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2. 用fortran编写peridynamics模型的求解程序，将计算得到的坐标信息，位移，速度，损伤等以vtk格式保存
peridynamics 相关的理论参考 Silling教授的文章以及其他相关的。
用无网格法对peridynamics 求解后，可以得到以下信息：
材料点的坐标信息：coord(i,1), coord(i,2), coord(i,3)，i从0到nodenumber(材料点总数)；1，2，3分别代表x，y，z
材料点标量信息：损伤 dam(i,1); i表示材料点编号，1没啥信息，标量；当然还有其他标量比如压力，温度等；
材料点的矢量信息：位移disp(i,1), disp(i,2), disp(i,3); 速度 vel(i,1), vel(i,2), vel(i,3); i同样是材料点的编号，1，2，3表示x，y，z
求解是按照时间步推进的，每一个时间下，都会得到以上所有信息，我们需要每一个时间步都记录下以上所有信息，每一个时间步对应一个vtk文档，
不同时间步对应的vtk文档的命名采用数字递增，这样在paraview中就可以识别到一个vtk群组，可形成动画效果。
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
fortran 代码：
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
implicit none
% 申明变量
integer:: i,j,k,tt,nt = 2000
character(len=100):: fileDir,fileName,fileSeq,fileSuffix,myStr
integer::counterLen = 4
logical::FEXIST
fileDir = 'D:\0923\' % fileDir 表示文件路径；这里是 d盘下面的0923文件夹；后续生成的.vtk文件都会保存在这里
inquire(file = trim(fileDir),EXIST=FEXIST) % 看看该文件夹是否已经存在
if(.not.FEXIST) then % 如果没有，就系统创建一个
call system('md '//trim(fileDir)) %
endif
Do k=tt, nt ! 注释 nt是总的时间步，例如2000步，根据对peridynamics的求解需要
fileName = 'view' % 文件名 为 view
fileSuffix = '.vtk' % 文件后缀为 .vtk
write(myStr,'(I2)')counterLen % 参数counterLen 是用来连续命名.vtk文件的，这里取 counterLen等于4，意味着最大将生成 9999个文件；若取为3，则只能生成999个；
%% 要保证每个时间步都对应一个.vtk文件，counterLen应该根据总时间步nt来确定； counterLen等于4，第一、二个文件是 view0001.vtk; view0002.vtk
myStr = '(I'//trim(adjustl(myStr))//')' % myStr决定文件名后面的是0001还是0002；
write(fileSeq,myStr)tt % 这样文件名后面就是tt了，也就是当前时间步；比如当前时间步为第100，也就是tt=100；那么对应的vtk文件就是 view0100.vtk
do j = 1,counterLen - len(trim(adjustl(fileSeq)))
fileSeq = '0'//adjustl(trim(fileSeq))
enddo
open(tt+50,file = adjustl(trim(fileDir))//adjustl(trim(fileName))&
&//adjustl(trim(fileSeq))//adjustl(trim(fileSuffix))) % 打开该生成的文件，准备写入了；
write(tt+50,168) %% 写入第一行，表明vtk的版本；第二行，标题，随便取；第三行，类型，我这里用的是 DATASET UNSTRUCTURED_grid 对应于只有材料点
write(tt+50,169) %% 写入 'POINTS',3x,'20000',3x,'float'
do i = 1, totnode
write(tt+50,170) coord(i,1), coord(i,2), coord(i,3) %% 写入点的坐标信息
end do
write(tt+50,171) %% 写入这句话 'point_data',3x,'20000',/,'SCALARS Damage FLOAT 1',/,'LOOKUP_TABLE DEFAULT' / 表示换行
do i = 1, totnode
write(tt+50,172) dmg(i,1) %% 写入 每个点的损伤信息
end do
write(tt+50,173) %% 写入 这句话 'VECTORS Displacement FLOAT'
do i = 1, totnode
write(tt+50,170) disp(i,1), disp(i,2), disp(i,3) %% 写入 点的位移矢量信息
end do
write(tt+50,174) %% 写入这句话'VECTORS Velocity FLOAT'
do i = 1, totnode
write(tt+50,170) vel(i,1), vel(i,2), vel(i,3) %% 写入每个点的速度信息
end do
close(tt+50)
end do
168 format('# vtk DataFile Version 3.0', / ,'title_impact', /, 'ASCII', /, 'DATASET UNSTRUCTURED_grid') %% 3x表示空三个位置；/表示换行；
169 format('POINTS',3x,'182709',3x,'float')
170 format(e12.5,3x,e12.5,3x,e12.5)
171 format('point_data',3x,'182709',/,'SCALARS Damage FLOAT 1',/,'LOOKUP_TABLE DEFAULT')
172 format(e12.5)
173 format('VECTORS Displacement FLOAT')
174 format('VECTORS Velocity FLOAT')
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
附生成的vtk格式文档压缩包
百度网盘链接：https://pan.baidu.com/s/1PjiJDZuRLH6i3u5gSkD9Cw
提取码：0104
整个有1350个文件，每个10M左右，太大了，因此我只放了最后的大概10个
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
3. 用paraview 看看生成的vtk格式文件
1. 用paraview打开生成的vtk文件；因为生成的VTK文件是连续编号的，paraview自动识别为组群。点击最上头的view...vtk 即可
2. 点击apply，把surface改为point gaussian ，因为我们的是点，没有cell；
3. 就可以看到点构成的结构啦。
4. 我这里主要想看损伤。点击左边选择 damage；点击edit，会出来右边，可以调节范围；在 rescale那里，选择 clamp and update every time step ;
点击右上角的小x，关掉右边的这一栏；
5. 这里可以播放，可以选择要显示的时间步。
6. 我们看看最后的效果吧
四、生成.VTK文件的子程序
! subroutine to output data to paraview readable vtk format
!
subroutine output_to_paraview_vtk(itimestep, x, vx, dis,dvx, mass, rho, p, u, c, itype, hsml, ntotal, mark,tots)
use parameter
implicit none
character(len=100):: fileDir,fileName,fileSeq,fileSuffix,myStr
integer::counterLen = 6
logical::FEXIST
integer:: i, j, itimestep,ntotal, itype(maxn), nvirt, mark(maxn)
double precision x(3, maxn), vx(3, maxn), dvx(3, maxn), dis(3, maxn),mass(maxn), rho(maxn), p(maxn), u(maxn), hsml(maxn), c(maxn),pforce(maxn,3),tots(maxn,3,3)
fileDir = 'D:\PD-SPH-paraview\'
! inquire(file = trim(fileDir),EXIST=FEXIST)
! if(.not.FEXIST) then
! call system('md '//trim(fileDir))
! end if
fileName = 'view'
fileSuffix = '.vtk'
write(myStr,'(I2)')counterLen
myStr = '(I'//trim(adjustl(myStr))//')'
write(fileSeq,myStr)itimestep
do j = 1,counterLen - len(trim(adjustl(fileSeq)))
fileSeq = '0'//adjustl(trim(fileSeq))
end do
open(itimestep+50,file = adjustl(trim(fileDir))//adjustl(trim(fileName))&
&//adjustl(trim(fileSeq))//adjustl(trim(fileSuffix)))
! vtk inofrmation
write(itimestep+50,168)
! particle position
write(itimestep+50,169)
do i = 1, ntotal
if(abs(itype(i)).ne.3) then
if (dim .eq. 1) then
write (itimestep+50, 170) x(1, i), 0.0, 0.0
else if (dim .eq. 2) then
write (itimestep+50, 170) x(1, i), x(2, i), 0.0
else if (dim .eq. 3) then
write (itimestep+50, 170) x(1, i), x(2, i), x(3, i)
end if
end if
end do
! particle mass
write(itimestep+50,171)
do i = 1, ntotal
if(abs(itype(i)).ne.3) then
write (itimestep+50, 172) mass(i)
end if
end do
! particle density
write(itimestep+50,173)
do i = 1, ntotal
if(abs(itype(i)).ne.3) then
write (itimestep+50, 172) rho(i)
end if
end do
! particle pressure
write(itimestep+50,174)
do i = 1, ntotal
if(abs(itype(i)).ne.3) then
write (itimestep+50, 172) p(i)
end if
end do
! particle particle_type
write(itimestep+50,178)
do i = 1, ntotal
if(abs(itype(i)).ne.3) then
write (itimestep+50, 180) itype(i)
end if
end do
! particle smooth_length
write(itimestep+50,179)
do i = 1, ntotal
if(abs(itype(i)).ne.3) then
write (itimestep+50, 172) hsml(i)
end if
end do
168 format('# vtk DataFile Version 3.0', / ,'title_impact', /, 'ASCII', /, 'DATASET UNSTRUCTURED_grid')
169 format('POINTS',3x,'584808',3x,'float')
170 format(e12.5,3x,e12.5,3x,e12.5)
171 format('point_data',3x,'584808',/,'SCALARS mass FLOAT 1',/,'LOOKUP_TABLE DEFAULT')
172 format(e12.5)
180 format(I5)
173 format('SCALARS density FLOAT 1',/,'LOOKUP_TABLE DEFAULT')
174 format('SCALARS pressure FLOAT 1',/,'LOOKUP_TABLE DEFAULT')
175 format('SCALARS energy FLOAT 1',/,'LOOKUP_TABLE DEFAULT')
176 format('VECTORS Velocity FLOAT')
177 format('VECTORS Acceleration FLOAT')
178 format('SCALARS particle_type FLOAT 1',/,'LOOKUP_TABLE DEFAULT')
179 format('SCALARS smooth_length FLOAT 1',/,'LOOKUP_TABLE DEFAULT')
close(itimestep+50)
end
