marc焊接仿真二次开发-以几何位置激活生死单元的方法？ 50
二次开发指南上和网上关于生死单元的二次开发uactive（）大多都是用单元号作的。这样在遍历单元格时很容易判断所遍历的单元生死状态，但对模型分网的网格单元号是有要求的。二次开发指南上网格单元号沿焊接路径方向有规律增加。可是对于不在marc上建模分网的焊接模型，网格号很难满足上述要求。
此时想用遍历单元与移动的热源中心相对位置来决定是否激活单元。我们看到二次开发指南中将生死单元状态另写入以文件中，仿照此方法，在初始化网格生死状态后，将所有会出现生死状态变化的单元的单元号和生死状态写入一.txt文件。进入焊接工况下，为判断遍历单元生死状态以决定生死操作，先将.txt文件读入二维向量，进行生死操作后删除原.txt文件，再将更新后的二维数组数据重新写入.txt文件。
这种方法效率极低，运行极慢。各位仿真大佬有没有效率更高的方法，十分感谢。另附现有代码。
另外，marc程序中是否有与fluent相似的数据结构，即用单元的结构体链表存储各个单元的各向参数。我看到elevar（）函数可以调出的项里没有单元的生死状态。
还有，弱弱问一句，死单元里有温度和热流等参数吗？
代码如下：
subroutine uactive(m,n,mode,irststr,irststn,inc,time,timinc)
include 'D:\MSC.Software\Marc\2016.0.0\marc2016\common\implicit'
dimension m(2)
real,dimension(3)::coordin !定义坐标
integer,dimension(35354,2)::modee
c
c 用来控制单元状态的子程序
c
c 输入参数:
c M(1) 单元号
c M(2) 在自适应分析中为母单元号
c N 积分点号
c INC 增量步号
c TIME 增量步开始的时间
c TIMINC 时间增量
c
c 输出参数:
c MODE: 单元的死活状态
c =-1,定义死单元
c =1,定义活单元
c =2,保持单元状态不变
c IRSTSTR:等于1时表示重新设置单元应力为0
c IRSTSTN:等于1时表示重新设置单元应变为0
c* * * * *
c
cc 添加三个公共变量,
cc foreInc保存增量步
cc iopen文件打开标志
cc iend计算完成标志
common/check/foreInc,iopen,iend
data iopen,iend/0,0/
real x0,y0,z0,x1,y1,z1,alpha,v,r,omega
c (x0,y0,z0)为圆心位置,(x1,y1,z1)为热源中心位置,v为激光扫描线速度，r为扫描路径半径，omega为激光扫描角速度，alpha为参数方程参量
c 在0-19.5s阶段，即第一层熔覆过程，设置圆心位置为修复层上表面圆心，
c 在19.5-43.5s阶段，即第二层熔覆过程，设置圆心位置为修复层上表面圆心
integer::ie,mode,icode,kcus,nn !ie 为单元号，icode=0直接返回坐标，连续体单元kcus=1
real::dist,dire !dist为遍历单元到热源中心距离，dire为数量积
ie=m(1);
if(ie.gt.25712.and.ie.lt.502035) return !set直接跳过
icode=0;kcus=1;nn=n
call elmvar(icode,m,nn,kcus,coordin)
cc
if(time.le.19.5) then
x0=175.54; y0=187.363;z0=-187.388
else if (time.gt.19.5.and.time.le.43.5)then
x0=178.72;y0=190.471; z0=-185.311
end if
c 设置激光扫描的线速度常数10,定义角速度omega
v=10;omega=v/r;pi=3.1415926
c 设置激光扫描热源中心位置参数方程重要参量alpha
if(time.le.13)then
r=19;alpha=omega*time
else if (time.gt.13.and.time.le.19.5)then
r=9;alpha=omega*time-13*omega
else if(time.gt.19.5.and.time.le.35.5)then
r=24;alpha=omega*time-19.5*time
else if (time.gt.35.5.and.time.le.43.5)then
r=12;alpha=omega*time-35.5*time
end if
c 设置激光热源中心位置，假定热源在法线（1，1，1）的平面上运动
x1=x0+(cos(alpha)-sin(alpha))*r/1.4142135
y1=y0+cos(alpha)*r/1.4142135
z1=z0+(cos(alpha)+sin(alpha))*r/1.73205081 热源位置
cc
dist=sqrt((coordin(1)-x1)**2+(coordin(2)-y1)**2
$+(coordin(3)-z1)**2) 积分点热源距离
dire=3.18*(175.54-coordin(1))+3.108*(187.363-coordin(2))
$+2.077*(-187.388-coordin(3)) 数量积确定位置
c 初始激活第一层热源周围0.5mm的单元 ,死单元存入for96
if (inc.eq.0) then
if(dist.le.0.5)then
if(dire.ge.0)then
mode=1
else mode=-1
end if
else mode=-1
end if
c 死单元初始化
open(95,file='mode.txt',position='append')
write(95,*)ie,mode
close(95)
go to 11
end if
open(95,file='mode.txt',status='old')
do im=1,35454
read(95,*) modee(im,1),modee(im,2)
if(ie.eq.modee(im,1)) then
mode=modee(im,2)
end if
end do
close (95)
if(time+timinc.le.19.5)then
if(mode.eq.-1)then
if(dist.le.11)then
if(dire.ge.0)then
mode=1
else mode=-1
end if
else mode=-1
end if
else
return
end if
end if
if(time+timinc.le.43.5.and.time+timinc.gt.19.5)then
if(mode.eq.-1)then
if(dist.le.12) then
mode=1
else mode=-1;
end if
else mode=-1;
end if
else
return
end if
c 一轮生死单元判断过后，更新mode.txt，并记录死单元
open(95,file='mode.txt',status='old')
close(95,status='delete')
c 删除原mode.txt文件，再重建更新
open(95,file='mode.txt',status='new')
do im=1,35454
write(95,*) modee(im,1),modee(im,2)
end do
close(95)
c
cc
c 将当前增量步的死单元号写入文件 fort.96 中
c 以便查看或者进行其他的处理
11 if(mode.eq.-1) then
write(96,*) 'deactivating element ', ie, ' increment ', inc
endif
cc
cc 新建一个过程文件 PostHide.proc 用来辅助后处理
if(iopen.eq.0) then
iopen=1
foreInc=inc
open(50,file='PostHide.proc',status='unknown')
close(50,status='delete')
open(50,file='PostHide.proc',status='new')
write(50,*)'*select_elements'
endif
if(iend.eq.0) then !判断是否到最后增量步
if(inc.eq.foreInc) then !判断增量步是否变化
if(mode.eq.-1) write(50,*) ie !写入当前增量步的死单元号
else
write(50,*)'#' !选择结束标志
write(50,*)'*invisible_selected' !使选中单元不可见
write(50,200) (inc-1) !显示第 (inc-1) 增量步
200 format('*post_skip_to',i5,/,'*select_elements')
write(50,*)'*animation_save' !保存当前窗口图形
cc
if(mode.eq.-1) write(50,*) ie
foreInc=inc
endif
if(inc.eq.435) then !到最后一步,置 iend=1
iend=1
write(50,*)'#'
write(50,*)'*invisible_selected'
write(50,*)'*animation_save'
endif
endif
cc
c* * * * *
return
end
代码尚未跑通，见笑了
