0032] (3)从存储邻接体单元的集合的哈希表中提取当前结点包含的体单元，将所有单 元访问位为true的单元压入当前单元分组，并设置这些单元的单元访问位为true，然后将 所包含单元其它结点中结点访问位为false的结点依次压入结点队列。
[0033] (4)判断所有节点是否完成遍历，满足条件转步骤6,否则判断单元数是否达到分 组单元数上限，条件满足则将当前分组设置为下一分组。
[0034] (5)从队列中提取第一个结点编号，并将该结点访问标识设置为true，转步骤3。
[0035] (6)分组完成后，将各组内结点根据所属单元列表是否完全在当前分组内进行分 组，满足条件的为组内点，否则为边界点。对所有结点依照分区及优先组内点，然后编号边 界点的顺序进行重编号。
[0036] (7)根据新的编号，更新所有的单元结点编号信息。并将结点坐标信息和单元结点 编号信息存储在一起。
[0037] (8)根据各单元顶点信息采用等参元法计算得到各单元和标准单元形函数的变换 关系矩阵，各结点质量矩阵及单元受力索引。其中等参元法是将各个四面体单元与正四面 体的形函数之间建立映射关系，表现形式为一个3*3矩阵，在每一个时间步中计算顶点坐 标发生变化后，通过该变换关系获取新的形函数矩阵。同时输入模型参数杨氏模量、泊松比 及时间步。
[0038] (10)GPU端初始化：将除结点坐标信息外的其它信息传送到GPU端。
[0039] (11)根据手动输入或是力反馈设备捕获的边界受力变化或是扫描仪或相机捕获 的可见表面节点的位移变化，利用Neo-Hookean模型，该模型定义了力与顶点位移之间的关 系，总体拉格朗日方法，计算每一个时间步的组成四面体单元的各顶点的受力。在多GPU的 并行实现中将四面体单元的顶点坐标信息分成若干段依次传输给GPU，同时在每段数据传
输完成后开启GPU端该段数据的计算函数 1 '，： 其中S表示第二Piola-Kirchoff应力，μ表示剪切模量，J表示变形梯度，I是三阶单位 矩阵，L是右柯西格林张量第一不变量，C是右柯西格林张量，δ是材料的体积弹性模量， 通过该公式计算各单元结点的受力情况。其中分段大小依据GPU性能动态设定。这种分段 处理方式可以减少多个GPU之间的数据交互。
[0040] (12)各单元结点受力计算完成后，通过预存的单元受力索引计算结点所受合力。
[0041] (13)将受力信息传送回CPU端之后，调用碰撞检测和响应程序进行处理后，采用 中心差分法公式1
各节点的坐标变化计算出下一 步的结点位移。其中At表示时间步，Μ表示四面体单元的质量，un+1表示下一个时间步的 位移，Rn表示当前时刻的反作用力，？"表示当前时刻所施加的外力，uni表示前一时刻的位 移，un表示当前时刻的位移即步骤二中捕获的表面节点的位移变化，更新四面体单元各顶 点的坐标变化。判断模型是否已经收敛稳定。条件满足则中止计算。否则转步骤11继续 计算模型变化。
[0042] 如图3所示，左侧为变形前的模型，右侧为变形后的结果。
[0043] 本发明中未详细阐述的部分属于本领域技术人员的公知技术。
【主权项】
1. 一种基于多GPU的人脑变形仿真方法，其特征在于包括如下步骤： (1) 初始化，读入人脑四面体网格模型数据，分析四面体网格单元的拓扑结构关系得到 四面体网格邻接体单元的集合，依据所获得的集合信息和所用GPU核数目，对读入的四面 体网格模型原始数据进行分组、编号，然后基于等参元法计算各四面体单元的形函数矩阵， 用于后面在发生受力或表面位移后内部单元的形变计算； (2) 初始化完成后，根据手动输入或是力反馈设备捕获的边界受力变化或是扫描仪或 相机捕获的可见表面节点的位移变化，利用Neo-Hookean模型，该模型定义了力与顶点位 移之间的关系，总体拉格朗日方法，计算每一个时间步的组成四面体单元的各顶点的受力， 在实现过程中采用多GPU加速计算过程，进行数据分组减少多GPU数据交互； (3) 把步骤二中每个GPU核计算出的受力结果进行合并，并根据最终计算结果采用中，其中At表示时间步，M表示四面 体单元的质量，un+1表示下一个时间步的位移，Rn表示当前时刻的反作用力，Fn表示当前时 刻所施加的外力即步骤二捕获的表面受力，Unl表示前一时刻的位移，Un表示当前时刻的位 移即步骤二中捕获的表面节点的位移变化，更新四面体单元各顶点的坐标变化，并在多CPU 上实现。2. 根据权利要求1所述的基于多GPU的人脑变形仿真方法，其特征在于：步骤⑴中 的对读入的四面体网格模型原始数据进行分组、编号，然后基于等参元法计算各四面体单 元的形函数矩阵，包括步骤如下： (2. 1)遍历所有四面体网格单元的顶点，为每一个顶点建立所属单元列表，为每个顶点 建立访问标志位并设置为false; (2. 2)输入分组起始参考顶点，对四面体网格模型单元进行分组，分组单元数=四面 体总单元数/所用GPU核数，首先将起始点的所属单元加入分组，并将顶点访问标志设为 true，判断当前分组内包含四面体是否达到分组单元数限制，达到则从下一个输入的四面 体单元顶点开始，将四面体放入下一分组，否则读取参考点所属单元列表内访问标志顶点 为false的节点的所属单元并加入当前分组中； (2.3) 步骤（2.2)对四面体单元进行分组完成后，将各组内顶点根据所属单元列表是 否完全在当前分组内进行区分，满足条件的为组内点，否则为边界点；对所有分组内四面体 单元的顶点依照分区次序从低到高，对这些四面体单元组内顶点和边界点进行重编号； (2.4) 根据新的编号，更新所有的单元顶点编号信息； (2. 5)根据各单元顶点信息采用等参元法计算得到各单元和标准单元形函数的变换关 系矩阵，等参元法是将各个四面体单元与正四面体的形函数之间建立映射关系，表现形式 为一个3*3矩阵，在每一个时间步中计算顶点坐标发生变化后，通过该变换关系获取新的 形函数矩阵。3. 根据权利要求1所述的多GPU的人脑变形仿真方法，其特征在于：步骤⑵中采用 多GPU加速计算过程，进行数据分组减少多GPU数据交互，包括步骤如下： (3. 1)初始化：初始化完成后，将四面体单元的顶点的坐标信息和单元节点编号信息 存储于同一数据结构中，保证在单元信息传送到GPU端后，已传送四面体单元相关的顶点 信息也已经传送完毕； (3. 2)流式传输：为了充分利用多GPU服务器上每个GPU核的计算能力及数据传输带 宽，把输入数据分成若干组，每次传输一组数据后，在开启GPU计算的同时进行下一组数据 的传输； (3. 3)在GPU端首先使用由Neo-Hookean模型及总体拉格朗日方法推导出的第二 Piola-Kirchoff应力其中S表示第二 Piola-KirchofT应力，μ表示剪切模量，J表示变形梯度，I是三阶单位矩阵，^是右柯西 格林张量第一不变量，C是右柯西格林张量，δ是材料的体积弹性模量，通过该公式计算各 单元的受力情况； (3. 4)通过初始化计算时生成的各四面体单元顶点受力求和求整个四面体网格模型各 顶点的总受力，并将结果传送回CPU端。4.根据权利要求1所述的多GPU的人脑变形仿真方法，其特征在于：所述步骤⑵中总 体拉格朗日方法实现过程为：方法在预测模型变形时通过把整个连续的变形过程离散为若 干个时刻的变形，设计算开始的时间为t，计算变形的下一个时刻和时间t的时间间隔叫做 时间步，记做At，首先使用顶点的初始坐标、步骤（1)的形函数矩阵及捕获的表面位移或 受力计算时刻t+At时每个四面体单元四个顶点的受力情况，然后通过步骤（3)的计算获 得t+Δt时刻顶点坐标、形函数矩阵及捕获的表面位移或受力计算时刻t+2*At时刻t+Δt 时每个四面体单元四个顶点的受力情况以此类推知道两个时间步之间顶点位移不再发生 变化。
【专利摘要】本发明涉及一种基于多GPU的人脑变形仿真方法，使用总体拉格朗日法、中心差分法及非线性有限元Neo-Hookean模型来进行脑部变形的显式迭代仿真，通过在GPU端的并行化算法实现加速仿真计算的速度，在引入多GPU进一步提高并行能力之后，对原始数据采用广度优先搜索方法进行数据分组和重新编号最小化各结点之间的数据相关性，并通过使用额外的数据结构使得单个节点内的输入数据及计算过程存在一一对应关系，进而可以使用流式传输实现计算和数据传输的并行化，最终利用混合CPU/GPU的多核计算架构进行加速，满足了人脑变形仿真中的精度与速度要求。本发明充分利用了人脑变形中的细微性以及时空连续性特点，可以在计算机构造的虚拟环境中逼真地模拟脑部的变形过程。
【IPC分类】G06T17/30, G06T1/20
【公开号】CN105389853
【申请号】CN201510731095
【发明人】胡勇, 田野, 沈旭昆
【申请人】北京航空航天大学
【公开日】2016年3月9日
【申请日】2015年11月2日