CN103970960A - 基于gpu并行加速的无网格伽辽金法结构拓扑优化方法 - Google Patents

Abstract

本发明公开了一种基于GPU并行加速的无网格伽辽金法结构拓扑优化方法。本发明主要步骤如下：(1)将数据读入主机内存，通过CPU布置积分点，建立节点、积分点与局部搜索小格子的关系，计算节点影响域半径和积分点定义域半径，确立节点和积分点关系，然后将数据复制到GPU全局存储器中；(2)根据不同的计算数据设定不同的GPU线程块与线程数量；(3)通过CPU和GPU异步组装和求解无网格伽辽金法总体离散系统方程，得到位移近似解；(4)在GPU中进行结构拓扑优化计算，通过设计变量的残差来判断迭代是否结束及结果是否输出。本发明硬件成本低，通用性强，在满足工程精度要求的前提下可降低大量时间消耗。

Description

基于GPU并行加速的无网格伽辽金法结构拓扑优化方法

技术领域

本发明涉及计算机辅助工程中的结构优化领域，具体涉及一种基于GPU(Graphics Processing Unit)并行加速的无网格伽辽金法结构拓扑优化方法。

背景技术

近二十年来，在计算力学领域，无网格法脱颖而出。无网格法不是采用对场变量插值和近似的网格单元，而是采用对场变量进行逼近和近似的节点，同时利用权函数来表征节点信息。现在已经发展了十余种无网格方法，如光滑质点流体动力学法、无网格伽辽金法等。无网格伽辽金法是一种非常流行的无网格法，它采用移动最小二乘法构造位移近似形函数，从能量泛函的弱变分形式得到整体控制方程，再利用拉格朗日法、罚函数法等处理本质边界条件，最后求解总体离散系统方程得到节点位移的数值解。无网格伽辽金法具有较高的协调性和稳定性，计算精度高，收敛快，而且消除了剪切锁死和体积锁死等缺点，现已被广泛应用于冲击动力学、金属塑性成型、结构优化等问题，但是其刚度矩阵的组装和总体离散系统方程求解复杂，且在计算节点影响域半径、组装刚度矩阵等时采用全局搜索，导致计算量大，速度慢，耗时长，这些缺点严重限制了无网格伽辽金法在结构拓扑优化问题中的应用。

结构拓扑优化在工程结构设计的初始阶段起着重要作用，可为设计者提供一个结构布局最优方案的概念性设计，现已被广泛应用于建筑、机械和航天等领域。结构拓扑优化是使设计域内的材料变为有或者无的过程，即对结构传递力作用贡献小的区域，应删除材料，而对结构传递力作用有益的部位，应保留材料。为求得最佳的材料分配方案，使材料能最大程度的利用，需要构造结构拓扑优化问题的优化模型，确立设计变量，约束条件，以及目标函数，经过计算求得目标函数的最小值。结构拓扑优化可以节约材料，降低成本，但是优化模型求解要经过多次迭代计算，耗时很长。

随着计算技术和计算机技术的迅猛发展，并行计算已经在多个领域内取得了成功，成为提高计算效率，节约计算成本的有效方法。并行计算是指在并行计算机上，将一个任务分解成多个子任务，然后分配给不同的处理器，各个处理器之间相互协同，并行地执行子任务，从而达到提高求解速度、增大求解问题规模的目的。现有的并行技术主要有分布式并行计算、并行机、多线程等，主要是以CPU(Central Processing Unit)为计算核心的计算平台，虽然可以取得较好的加速比，但同时存在以下的不足，首先，以CPU为计算核心的并行机价格昂贵，维护和使用相对复杂，对算法的设计和编程能力要求也较高，较难普及；其次，对于大规模问题，会导致分布式计算机之间或者多线程之间的通讯消耗提升，从而降低了整体加速比。

2006年，NVIDIA公司公布了统一着色器模型和GeForce8系列GPU，GPU从此进入了通用计算机时代。统一着色器模型整合了顶点着色器和片段着色器，使每个着色器都可以自由承担原本某种特定的着色器的工作，这种方式使得GPU在图形处理时空闲的着色器更少，计算效率更高。随后，NVIDIA公布了统一计算设备架构CUDA(Compute Unified Device Architecture)，并发布了专门针对CUDA架构的编程语言CUDA C。CUDA C支持C语言的大多数指令和语法，并加入了使程序能在GPU上进行多线程计算的扩展语言。NVIDA公司推出的CUDA平台，作为一种通用的GPU并行计算环境，已成为当今最流行的GPU并行计算平台。

随着GPU的广泛应用，其强大的浮点运算能力越来越受到人们的重视，基于GPU的通用并行计算技术得到快速发展。在数值代数中，存在较为密集的数据计算，如矢量点积、线性方程组求解等，当处理的规模较大时，计算量是非常可观的，并且大多数科学与工程计算问题，最后的步骤一般归结为求解线性方程组，因此研究GPU加速代数计算，尤其是大规模线性方程组的GPU加速求解，其加速比一般可达到15～25倍，但是目前无网格伽辽金法结构拓扑优化仍然不能利用基于CUDA架构的GPU通用计算平台来实现并行加速。

发明内容

本发明的目的在于针对无网格伽辽金法的结构拓扑优化中存在计算量大、耗时长等问题，提供一种基于GPU并行加速的无网格伽辽金法结构拓扑优化方法。

本发明基于GPU并行加速的无网格伽辽金法结构拓扑优化方法，包括如下顺序的步骤：

(1)将无网格计算模型的节点坐标、有限元背景积分网格数据、局部搜索格子数据、本质边界数据、惩罚系数、节点影响域系数、材料物理参数、结构拓扑优化参数导入主机内存中；

(2)在CPU中根据主机内存中存储的模型有限元背景积分网格数据、本质边界数据，分别为模型内部、力边界和位移边界布置积分点，当为二维几何模型问题时，利用程序自带的Hammer和Gauss数值积分信息，在模型的有限元三角形背景积分网格和线边界上布置积分点，当为三维几何模型问题时，利用程序自带的Hammer数值积分信息，在模型的有限元四面体背景积分网格和本质边界上的有限元三角形背景积分网格中布置积分点，求得到积分点的坐标、权系数和雅可比系数；

(3)根据局部搜索格子数据划分出每个局部搜索小格子的区域，然后分别求得模型节点、积分点和局部搜索小格子的关系数据；

(4)采用局部搜索法计算每个节点的影响域半径和每个积分点的定义域半径，然后计算每个节点影响域内的积分点数据，以及每个积分点定义域内的节点数据；

(5)将模型节点坐标、模型内部和位移边界上的积分点数据、惩罚系数、材料物理参数、节点与局部搜索小格子的关系数据、积分点与局部搜索小格子的关系数据、结构拓扑优化参数、节点影响域内的积分点数据、积分点定义域内的节点数据复制到GPU全局存储器中；

(6)通过GPU读取GPU全局存储器中存储的模型节点坐标、积分点数据、积分点与局部搜索小格子的关系数据、积分点定义域内的节点数据，然后设置GPU中的线程块的数量和每一个线程块所包含的线程数量，并行加速计算每个积分点定义域内各个节点相应的形函数值及形函数导数值；

(7)根据GPU全局存储器中存储的积分点数据、积分点定义域内的节点数据、节点影响域内的积分点数据、结构拓扑优化参数、节点形函数值，首先定义模型节点的初始相对密度参数，然后计算结构总体积的初始值及灵敏度；

(8)在GPU执行步骤(6)～(7)时，根据模型节点数据、节点影响域半径、节点与局部搜索小格子的关系数据、局部搜索格子数据，在CPU中确定影响域有重合区域的两两节点，组成交叉节点对，并将交叉节点对信息存储于数组中；

(9)在CPU中读取主机内存中存储的模型节点坐标、本质边界上的积分点数据、本质边界数据、惩罚系数，利用罚函数法施加本质边界条件，在CPU中求得总体外力向量和总体惩罚力向量，然后将交叉节点对数据、总体外力向量、总体惩罚力向量复制到GPU全局存储器中；

(10)通过GPU进行无网格伽辽金法结构拓扑优化的迭代计算；

(11)判断是否需要输出无网格伽辽金法结构拓扑优化的计算结果，如果需要则将GPU全局存储器中存储的结构拓扑优化的计算结果复制到主机内存中，再写入输出文件。

具体地说，所述步骤(10)包括如下步骤：

(a)通过GPU读取GPU全局存储器中存储的节点形函数值、模型内部积分点数据、积分点定义域内的节点关系数据、模型内部积分点处的初始相对密度参数，然后建立线程与积分点一一映射模式，并行加速计算模型内部积分点处的相对密度值；

(b)通过GPU读取GPU全局存储器中存储的模型节点坐标、节点的影响域半径、材料物理参数、积分点数据、节点形函数值及形函数导数值、节点与局部搜索小格子的关系数据、积分点与局部搜索小格子的关系数据、本质边界数据、惩罚系数，然后建立GPU线程块与交叉节点对一一映射模式，GPU线程和组成交叉节点对的两个节点公共影响域内的积分点一一映射模式，并行加速组装总体刚度矩阵和总体惩罚刚度矩阵，并以按行压缩存储CSR格式存储于GPU全局存储器中；

(c)通过GPU读取GPU全局存储器中存储的总体刚度矩阵、总体惩罚刚度矩阵、总体外力向量和总体惩罚力向量，调用CUDA架构平台提供的免费函数库CUBLAS，采用共轭梯度法求解总体离散系统方程，求得无网格伽辽金法的位移近似解；

(d)通过GPU读取GPU全局存储器中存储的积分点数据、节点影响域内的积分点数据、积分点定义域内的节点数据、节点形函数值及形函数导数值、节点无网格伽辽金法位移近似解、材料物理参数以及结构拓扑优化参数，建立线程块与节点一一映射模式，线程与节点影响域内的积分点一一映射模式，并行计算目标函数的灵敏度；

(e)利用OC准则，通过GPU读取GPU全局存储器中的模型内部积分点处的相对密度值、结构总体积的初始值和灵敏度、目标函数的灵敏度、节点形函数值、积分点数据、积分点定义域内的节点数据，求得相对密度参数和体积比；

(f)计算设计变量前后两次迭代的残差，并判断残差是否满足初始残差限制，如果满足则迭代结束，如果不满足，则重新开始步骤(a)。

本发明与现有技术相比，具有下述优点：

(1)本发明采用三角形积分网格(二维几何问题)和四面体积分网格(三维几何问题)作为无网格伽辽金法的背景积分网格，三角形积分网格和四面体积分网格具有强大的适应性，使得本发明方法可广泛适用于二维、三维任意形状的几何模型问题。

(2)本发明利用局部搜索方法，分别建立模型节点、积分点与局部搜索小格子的联系，在进行无网格伽辽金法结构拓扑优化计算时，避免了费时的全局搜索，提高了计算效率。

(3)本发明将无网格伽辽金法总体刚度矩阵的组装和本质边界条件的施加结合起来，减少了总体惩罚刚度矩阵的组装时间，提高了计算效率，同时简化了程序，增加了程序的易读性和可移植性，另外本发明将总体刚度矩阵和总体惩罚刚度矩阵以按行压缩存储CSR格式存储于GPU全局存储器中，大大降低了对GPU全局存储器存储容量的需求。

(4)本发明将无网格伽辽金法模型节点形函数值及其形函数导数值的计算，总体离散方程的组装与求解，结构拓扑优化计算等均移植到GPU上进行并行加速计算，而积分点的布置，总体外力向量和总体惩罚力向量计算，数据的读入和输出等均在CPU上执行，而GPU和CPU采用异步执行的方式，有效的避免了由于主机内存和GPU全局内存之间数据交换所产生的额外时间。

(5)本发明在采用GPU进行计算时，根据并行计算对像和计算任务的不同，建立不同的线程映射模式，即主要建立三种映射模式：线程与节点一一映射，线程与积分点一一映射，线程块与交叉节点对一一映射。

(6)本发明能够充分利用GPU强大的浮点计算能力，为保证程序的计算效率、简洁性、可移植性，采用CUDA C编程，同时能够利用NVIDIA公司提供的免费线性代数数值计算库CUBLAS，所以本方法计算成本低，通用性强，程序简洁，具有很好的适应性。

(7)本发明通过建立GPU线程、线程块的各种映射模式，在保证求解精度的前提下，能够大幅度节约计算时间，提高计算效率，增大求解规模。利用基于CUDA构架的GPU强大的浮点计算能力，大规模的GPU计算线程，在进行无网格伽辽金法结构拓扑优化计算时一个线程块或一个线程负责一对交叉节点对、一个节点、一个积分点的相应计算，充分利用了具有上千处理核的GPU并行计算能力。

附图说明

图1为本发明方法的主流程框图。

图2为本发明方法中计算模型节点与局部搜索小格子的关系数据的流程框图。

图3为本发明方法中计算无网格伽辽金法结构拓扑优化目标函数的灵敏度的流程框图。

图4为本发明实施例的模型示意图

图5为图4所示实施例的模型节点示意图。

图6为图4所示实施例的背景积分网格示意图。

图7为图4所示实施例的结构拓扑优化结果示意图。

具体实施方式

下面结合附图和实施例对本发明作进一步详细的描述。

参见图1，本发明基于GPU并行加速的无网格伽辽金法结构拓扑优化方法的实施步骤如下：

(1)将无网格计算模型的节点坐标、有限元背景积分网格数据、局部搜索格子数据、本质边界数据、惩罚系数、节点影响域系数、材料物理参数、结构拓扑优化参数导入主机内存中；

(2)在CPU中根据主机内存中存储的模型有限元背景积分网格数据、本质边界数据，分别为模型内部、力边界和位移边界布置积分点，当为二维几何模型问题时，利用程序自带的Hammer和Gauss数值积分信息，分别在模型的有限元三角形背景积分网格和线边界上布置积分点，当为三维几何模型问题时，利用程序自带的Hammer数值积分信息，在模型的有限元四面体背景积分网格和本质边界上的有限元三角形背景积分网格中布置积分点，求得到积分点的坐标、权系数及雅可比系数；

(3)参见图2，根据局部搜索格子数据划分出每个局部搜索小格子的区域，然后分别求得模型节点、积分点和局部搜索小格子的关系数据；具体的步骤如下：

(a)根据局部搜索格子数据，划分局部搜索总格子，得到每个局部搜索小格子的范围，并以数据形式为每个小格子编号；

(b)根据模型节点编号从小到大顺序逐个读取主机内存中存储的模型节点坐标；

(c)根据模型节点的坐标和局部搜索小格子的范围，确定该节点所在小格子的编号；

(d)为cellnumnod数组中与该小格子的编号相对应处的元素值加1；

(e)判断模型节点是否编历，如果不是则跳到步骤(b)继续提取，如果是则继续下一步计算；

(f)根据局部搜索小格子的编号和cellnumnod数组中的元素值，计算每个局部搜索小格子内的第一个节点在cellnode数组中的下标值并存储于cellnode_idx数组中；

(g)判断是否编历局部搜索小格子，如果没有则跳至步骤(f)继续执行，如果已经编历，则执行下一步；

(h)根据模型节点编号从小到大顺序依次提取节点的坐标和局部搜索小格子的范围数据，确定该节点所在小格子的编号，并将该节点编号存储于数据cellnode中的相应位置；

(i)判断是否编历模型节点，如果没有则继续执行步骤(h)，如果编历则继续下一步。

(4)采用局部搜索法计算每个节点的影响域半径和每个积分点的定义域半径，然后确定每个节点影响域内的积分点数据，以及每个积分点定义域内的节点数据；

(5)将模型节点坐标、模型内部和位移边界上的积分点数据、惩罚系数、材料物理参数、节点与局部搜索小格子的关系数据、积分点与局部搜索小格子的关系数据、结构拓扑优化参数、节点影响域内的积分点数据、积分点定义域内的节点数据复制到GPU全局存储器中；

(6)通过GPU读取GPU全局存储器中存储的模型节点坐标、积分点数据、积分点与局部搜索小格子的关系数据、积分点定义域内的节点数据，设置GPU中的线程块的数量和每一个线程块所包含的线程数量，并行加速计算每个积分点定义域内各个节点相应的形函数值及形函数导数值；

(7)根据GPU全局存储器中的积分点数据、积分点定义域内的节点数据、节点影响域内的积分点数据、结构拓扑优化参数、节点形函数值，首先定义模型节点的初始相对密度参数，然后计算结构总体积的初始值和灵敏度；

(8)在GPU执行步骤(6)～(7)时，根据模型节点坐标、节点的影响域半径、节点与局部搜索小格子的关系数据、局部搜索格子数据，在CPU中确定影响域有重合区域的两两结点，组成交叉节点对，并将交叉节点对信息存储于数组中；

(9)在CPU中读取主机内存中存储的模型节点坐标、本质边界积分点数据、本质边界数据、惩罚系数，利用罚函数法施加本质边界条件，在CPU中计算总体外力向量和总体惩罚力向量，然后将交叉节点对数据、总体外力向量、总体惩罚力向量复制到GPU全局存储器中；

(10)通过GPU进行无网格伽辽金法结构拓扑优化的迭代计算；

(11)判断是否需要输出无网格伽辽金法结构拓扑优化的计算结果，如果需要则将GPU全局存储器上的结构拓扑优化结果复制到主机内存中，然后通过主机内存写入输出文件。

所述步骤(10)的详细步骤包括：

(a)通过GPU读取GPU全局存储器中存储的节点形函数值、模型内部积分点数据、积分点定义域内的节点数据、模型内部积分点处的初始相对密度参数值，建立线程与积分点一一映射模式，并行加速计算模型内部积分点处的相对密度值；

(b)通过GPU读取GPU全局存储器中存储的模型节点坐标、节点的影响域半径、材料物理参数、积分点数据、节点形函数值及形函数导数值、节点与局部搜索小格子的关系数据、积分点与局部搜索小格子的关系数据、本质边界数据、惩罚系数，然后建立GPU线程块与交叉节点对一一映射模式，GPU线程和组成交叉节点对的两个节点公共影响域内的积分点一一映射模式，加速组装总体刚度矩阵和总体惩罚刚度矩阵，并以按行压缩存储CSR格式存储于GPU全局存储器中；

(c)通过GPU读取GPU全局存储器中存储的总体刚度矩阵、总体惩罚刚度矩阵、总体外力向量、总体惩罚力向量，调用CUDA平台提供的免费函数库GUBLAS，采用共轭梯度法求解总体离散系统方程，得到模型各个节点的无网格伽辽金法位移近似解；

(d)参见图3，通过GPU读取GPU全局存储器中存储的积分点数据、节点影响域内的积分点数据、积分点定义域内的节点数据、节点形函数值及形函数导数值、节点无网格伽辽金法的位移近似解、材料物理参数以及结构拓扑优化参数，建立线程块与节点一一映射模式，线程与节点影响域内的积分点一一映射模式，并行加速计算目标函数的灵敏度；具体步骤如下：

第一步：通过GPU从按小到大逐个提取GPU全局存储器中存储的模型节点编号，根据节点编号和节点影响域内积分点数据得到该节点影响域内的积分点数和积分点编号；

第二步：依次提取该节点影响域内的积分点，根据积分点编号和积分点定义域内的节点数据，计算得到该节点影响域内的积分的定义域内的节点编号；

第三步：通过GPU读取GPU全局存储器中存储的节点形函数值、惩罚因子、内部积分点处的相对密度值、材料物理参数、积分点的权系数和雅可比系数、积分点的编号，计算该积分点处相对密度的(惩罚因子-1)次幂；

第四步：编历该节点影响域内的积分点定义域内的节点，根据该节点的编号和GPU全局存储器中存储的节点形函数导数值、节点位移近似值、该积分点处相对密度的(惩罚因子-1)次幂值，计算灵敏度分量值；

第五步：判断是否编历该节点影响域内的积分点定义域内的节点，如果没有，则跳至第四步继续执行，如果编历，则执行下一步；

第六步：累加灵敏度分量，得到灵敏度值；

第七步：判断是否编历该模型节点内的积分点，如果没有，则跳至第二步继续执行，如果编历则执行下一步；

第八步：判断是否编历模型节点，如果没有则跳至第一步继续执行，如果编历则执行下一步。

(e)利用OC(Optimization Criteria)准则，通过GPU读取GPU全局存储器中存储的模型内部积分点处的相对密度值、结构总体积的初始值和灵敏度、目标函数的灵敏度、节点形函数值、积分点数据、积分点定义域内的节点数据，求得相对密度参数和体积比；

(f)计算设计变量的前后两次迭代的残差，并判断残差是否满足初始残差限制，如果满足则迭代结束，如果不满足，则重新开始步骤(a)。

下面是将本发明方法应用于工程上的一个实例：

参见图4，本实施例为长15m，宽10m，厚度为1m，受集中力P＝1N作用的悬臂梁，材料的弹性E＝1.0，泊松比μ＝0.3，无网格伽辽金法结构拓扑优化模型包括1033个节点，其分布如图5所示，有限元背景积分网格包括1033个节点和1944个三角形网格，如图6所示。本发明方法针对本实施例的具体实施步骤如下：

第一步：将已经准备好的模型数据(包括模型的节点坐标、有限元背景积分网格顶点的编号和坐标、本质边界节点对中两节点的编号、本质边界信息、惩罚系数)、材料物理参数(包括材料弹性模量、泊松比)、节点影响域系数、局部搜索格子数据、拓扑优化参数(惩罚因子、最大体积保留率)输入到input_para.txt文件中，在程序运行时，系统自动将以上数据读入到主机内存中。一个标准的input_para.txt文件格式如下表一：

表一

第二步：依据主机内存中存储的局部搜索总格子边长和x轴方向的等分数，计算每个局部搜索小格子的边长；本实施例中局部搜索总格子x轴方向的边长为15m，y轴方向的边长为10m，x轴方向的等分数为30，所以本实施例中每个小格子x轴方向的边长为0.5m，y轴方向的边长可根据局部搜索总格子边长的等分数与边长的长度成比例求得约为0.33333m。

第三步：在CPU中根据本实施例三角形有限元背景积分网格顶点的坐标与编号、程序自带的Hammer数值积分数据，在模型内部布置积分点，求得积分点的坐标、权系数及雅可比系数，用于存储模型内部积分点信息的数组需要的存储空间大小为：1944(本实施例中有限元背景积分网格个数)*3(每个三角形有限元背景积分网格内布置的积分点数)*4*sizeof(double)字节。

第四步：通过CPU读取主机内存中存储的有限元背景积分网格的顶点坐标与编号，本质边界节点对中两节点的编号，根据程序自带的Gauss积分数据，在本质边界上布置积分点，求得位移和力边界上积分点的坐标、权系数、雅克比系数，现在分配给位移边界积分点数组的存储空间为：24(本实施例中位移边界节点对个数)*4(每个单元节点对上布置的积分点数)*4*sizeof(double)字节，因为本实施例中施加的是集中力，只在受力点布置一个积分点，所以为力边界积分点数组分配的存储空间大小为：4*sizeof(double)字节。

第五步：如图2所示，在CPU中，根据模型节点坐标、局部搜索格子数据，建立节点和局部搜索小格子的对应关系，并将其关系数据存储于cellnode，cellnumnod，cellnode_idx三个数组中，同理建立模型积分点与局部搜索小格子的对应关系。

第六步：在CPU中，对模型节点进行循环求得节点的动态影响域半径。根据节点坐标，确定该节点邻近的几个局部搜索小格子在全部格子中的编号，并根据局部搜索小格子与节点的对应关系数据确定该节点几个邻近局部搜索小格子中所有节点的编号，再在几个邻近小格子中进行全局搜索，求得邻近局部搜索小格子中其它节点和该节点的最小距离r_min，然后利用公式：节点影响域半径＝3.0(节点影响域系数)*r_min可求得该节点的动态影响域半径，在主机内存中为其分配的存储空间为：1033(模型节点的个数)*2(x、y轴两个方向)*sizeof(double)字节，同理求得积分点的定义域半径。

第七步：在CPU中，读取模型的节点坐标、积分点坐标、节点与局部搜索小格子的关系数据、积分点与局部搜索小格子的关系数据、节点的影响域半径，采用局部搜索方法求得每个节点影响域内的积分点数据，每个积分点定义域内的节点数据。

第八步：将模型节点坐标、模型内部和位移边界的积分点数据、位移边界数据、节点与局部搜索小格子的关系数据、积分点与局部搜索小格子的关系数据、结构拓扑优化参数、节点影响域内的积分点数据、积分点定义域内的节点数据复制到GPU全局存储器中，并为它们分配和主机内存中同样大小的存储空间。

第九步：通过GPU读取GPU全局存储器中存储的模型节点坐标、积分点数据、积分点与局部搜索小格子的关系数据、积分点定义域内的节点数据，建立线程与积分点一一映射模式，通过最小二乘逼近，并行求得每个积分点定义域内节点的形函数值，以及形函数对应于各个坐标轴的导数值，在本实施例中为形函数对x、y轴的导数值。

第十步：根据GPU全局存储器中存储的积分点数据、积分点定义域内的节点数据、节点影响域内的积分点数据、结构拓扑优化参数、节点形函数值，首先建立线程和节点一一映射模式，定义节点的初始相对密度参数，然后建立线程和积分点一一映射模式，算得结构总体积的初始值和灵敏度。

第十一步：GPU执行步骤九、十时，在CPU中根据节点坐标、节点的影响域半径、节点与局部搜索小格子的关系数据、局部搜索格子数据，采用局部搜索法，在CPU中确定影响域有重合区域的两两节点，组成交叉节点对，然后将交叉节点对信息以数据形式存储于主机内存中。

第十二步：读取主机内存中存储的模型节点坐标、本质边界数据、有限元背景积分网格数据、力边界积分点信息、节点与局部搜索小格子的关系数据、积分点与局部搜索小格子的关系数据，采用局部搜索法，求得模型的总体外力向量F。

第十三步：在CPU中读取主机内存中存储的模型节点坐标、位移边界积分点数据、本质边界数据、惩罚系数、节点与局部搜索小格子的关系数据、积分点与局部搜索小格子的关系数据，利用罚函数法施加本质边界条件，求得总体惩罚力向量F^a，然后将交叉节点对信息、总体外力向量和总体惩罚力向量复制到GPU全局存储器中。

第十四步：读取GPU全局存储器中模型节点的形函数值及形函数导数值、积分点数据、材料相对密度值，利用最小二乘逼近，在GPU中建立线程与积分点一一映射模式，计算模型内部积分点处的相对密度值。

第十五步：在GPU中建立线程块和交叉节点对一一映射模式，线程和组成交叉节点对的两个节点公共影响域内的积分点一一映射模式，并行加速组装总体刚度矩阵K和总体惩罚刚度矩阵K^a，并以CSR格式存储于GPU全局存储器中，然后读取GPU全局存储器中存储的总体外力向量F和总体惩罚力向量F^a，在GPU中调用CUDA平台提供的免费函数，采用共轭梯度法求解控制方程(K+K^a)U＝F+F^a，求得无网格伽辽金法的位移近似解U。

第十六步：如图3所示，在GPU中计算目标函数的灵敏度值。首先建立线程块与节点一一映射模式，求得每个节点影响域内的积分点数以及积分点编号，然后建立线程块内的线程与该节点影响域内的积分点一一映射模式，求得目标函数的灵敏度值。

第十七步：读取GPU全局存储器中存储的模型节点初始相对密度参数、结构总体积的初始值和灵敏度值、目标函数的灵敏度值、节点的形函数值、积分点数据、积分点定义域内的节点数据，建立积分点与GPU线程一一映射模式，在GPU上利用OC准则计算节点的相对密度值。

第十八步：在GPU上求得模型节点相对密度的残差，并判断残差是否小于初始设定值，如果大于则跳到第十四步继续迭代，如果小于则跳出迭代计算，并将GPU全局存储器中存储的无网格伽辽金法结构拓扑优化的计算结果复制到主机内存中，再写入输出文件，本实施例无网格伽辽金法结构拓扑优化结果如图7所示。

虽然参考优先实施例对本发明进行描述，但以上所述并不限定本发明的保护范围，凡属于本发明思路下的技术方案均属于本发明的保护范围，任何在本发明的精神及原则内的修改、改进等，也应视为本发明的保护范围。

Claims (4)

1.一种基于GPU并行加速的无网格伽辽金法结构拓扑优化方法，其特征在于包括如下顺序的步骤：

(1)将无网格计算模型的节点坐标、有限元背景积分网格数据、局部搜索格子数据、本质边界数据、惩罚系数、节点影响域系数、材料物理参数、结构拓扑优化参数导入主机内存中；

(2)在CPU中根据主机内存中存储的模型有限元背景积分网格数据、本质边界数据，分别为模型内部、力边界和位移边界布置积分点，当为二维几何模型问题时，利用程序自带的Hammer和Gauss数值积分信息，在模型的有限元三角形背景积分网格和线边界上布置积分点，当为三维几何模型问题时，利用程序自带的Hammer数值积分信息，在模型的有限元四面体背景积分网格和本质边界上的有限元三角形背景积分网格中布置积分点，求得到积分点的坐标、权系数和雅可比系数；

(3)根据局部搜索格子数据划分出每个局部搜索小格子的区域，然后分别求得模型节点、积分点和局部搜索小格子的关系数据；

(4)采用局部搜索法计算每个节点的影响域半径和每个积分点的定义域半径，然后计算每个节点影响域内的积分点数据，以及每个积分点定义域内的节点数据；

(5)将模型节点坐标、模型内部和位移边界上的积分点数据、惩罚系数、材料物理参数、节点与局部搜索小格子的关系数据、积分点与局部搜索小格子的关系数据、结构拓扑优化参数、节点影响域内的积分点数据、积分点定义域内的节点数据复制到GPU全局存储器中；

(6)通过GPU读取GPU全局存储器中存储的模型节点坐标、积分点数据、积分点与局部搜索小格子的关系数据、积分点定义域内的节点数据，然后设置GPU中的线程块的数量和每一个线程块所包含的线程数量，并行加速计算每个积分点定义域内各个节点相应的形函数值及形函数导数值；

(7)根据GPU全局存储器中存储的积分点数据、积分点定义域内的节点数据、节点影响域内的积分点数据、结构拓扑优化参数、节点形函数值，首先定义模型节点的初始相对密度参数，然后计算结构总体积的初始值及灵敏度；

(8)在GPU执行步骤(6)～(7)时，根据主机内存中存储的模型节点坐标、节点影响域半径、节点与局部搜索小格子的关系数据、局部搜索格子数据，在CPU中确定影响域有重合区域的两两节点，组成交叉节点对，并将交叉节点对信息存储于数组中；

(9)在CPU中读取主机内存中存储的模型节点坐标、本质边界上的积分点数据、本质边界数据、惩罚系数，利用罚函数法施加本质边界条件，在CPU中求得总体外力向量和总体惩罚力向量，然后将交叉节点对数据、总体外力向量、总体惩罚力向量复制到GPU全局存储器中；

(10)通过GPU进行无网格伽辽金法结构拓扑优化的迭代计算；

(11)判断是否需要输出无网格伽辽金法结构拓扑优化的计算结果，如果需要则将GPU全局存储器中存储的结构拓扑优化的计算结果复制到主机内存中，再写入输出文件。

2.根据权利要求1所述基于GPU并行加速的无网格伽辽金法结构拓扑优化方法，其特征在于所述步骤(10)包括如下步骤：

(a)通过GPU读取GPU全局存储器中存储的节点形函数值、模型内部积分点数据、积分点定义域内的节点关系数据、模型内部积分点处的初始相对密度参数，然后建立线程与积分点一一映射模式，并行加速计算模型内部积分点处的相对密度值；

(b)通过GPU读取GPU全局存储器中存储的模型节点坐标、节点的影响域半径、材料物理参数、积分点数据、节点形函数值及形函数导数值、节点与局部搜索小格子的关系数据、积分点与局部搜索小格子的关系数据、本质边界数据、惩罚系数，然后建立GPU线程块与交叉节点对一一映射模式，GPU线程和组成交叉节点对的两个节点公共影响域内的积分点一一映射模式，并行加速组装总体刚度矩阵和总体惩罚刚度矩阵，并以按行压缩存储CSR格式存储于GPU全局存储器中；

(c)通过GPU读取GPU全局存储器中存储的总体刚度矩阵、总体惩罚刚度矩阵、总体外力向量和总体惩罚力向量，调用CUDA架构平台提供的免费函数库CUBLAS，采用共轭梯度法求解总体离散系统方程，求得无网格伽辽金法的位移近似解；

(d)通过GPU读取GPU全局存储器中存储的积分点数据、节点影响域内的积分点数据、积分点定义域内的节点数据、节点形函数值及形函数导数值、节点无网格伽辽金法位移近似解、材料物理参数以及结构拓扑优化参数，建立线程块与节点一一映射模式，线程与节点影响域内的积分点一一映射模式，并行计算目标函数的灵敏度；

(e)利用OC准则，通过GPU读取GPU全局存储器中的模型内部积分点处的相对密度值、结构总体积的初始值和灵敏度、目标函数的灵敏度、节点形函数值、积分点数据、积分点定义域内的节点数据，求得相对密度参数和体积比；

(f)计算设计变量前后两次迭代的残差，并判断残差是否满足初始残差限制，如果满足则迭代结束，如果不满足，则重新开始步骤(a)。

3.根据权利要求2所述基于GPU并行加速的无网格伽辽金法结构拓扑优化方法，其特征在于所述步骤(d)包括如下步骤：

第一步：通过GPU从按小到大逐个提取GPU全局存储器中存储的模型节点编号，根据节点编号和节点影响域内积分点数据得到该节点影响域内的积分点数和积分点编号；

第二步：依次提取该节点影响域内的积分点，根据积分点编号和积分点定义域内的节点数据，计算得到该节点影响域内的积分的定义域内的节点编号；

第三步：通过GPU读取GPU全局存储器中存储的节点形函数值、惩罚因子、内部积分点处的相对密度值、材料物理参数、积分点的权系数和雅可比系数、积分点的编号，计算该积分点处相对密度的(惩罚因子-1)次幂；

第四步：编历该节点影响域内的积分点定义域内的节点，根据该节点的编号和GPU全局存储器中存储的节点形函数导数值、节点位移近似值、该积分点处相对密度的(惩罚因子-1)次幂值，计算灵敏度分量值；

第五步：判断是否编历该节点影响域内的积分点定义域内的节点，如果没有，则跳至第四步继续执行，如果编历，则执行下一步；

第六步：累加灵敏度分量，得到灵敏度值；

第七步：判断是否编历该模型节点内的积分点，如果没有，则跳至第二步继续执行，如果编历则执行下一步；

第八步：判断是否编历模型节点，如果没有则跳至第一步继续执行，如果编历则执行下一步。

4.根据权利要求1所述基于GPU并行加速的无网格伽辽金法结构拓扑优化方法，其特征在于所述步骤(3)包括如下步骤：

(a)根据局部搜索格子数据，划分局部搜索总格子，得到每个局部搜索小格子的范围，并以数据形式为每个小格子编号；

(b)根据模型节点编号从小到大顺序逐个读取主机内存中存储的模型节点坐标；

(c)根据模型节点的坐标和局部搜索小格子的范围，确定该节点所在小格子的编号；

(d)为cellnumnod数组中与该小格子的编号相对应处的元素值加1；

(e)判断模型节点是否编历，如果不是则跳到步骤(b)继续提取，如果是则继续下一步计算；

(f)根据局部搜索小格子的编号和cellnumnod数组中的元素值，计算每个局部搜索小格子内的第一个节点在cellnode数组中的下标值并存储于cellnode_idx数组中；

(g)判断是否编历局部搜索小格子，如果没有则跳至步骤(f)继续执行，如果已经编历，则执行下一步；

(h)根据模型节点编号从小到大顺序依次提取节点的坐标和局部搜索小格子的范围数据，确定该节点所在小格子的编号，并将该节点编号存储于数据cellnode中的相应位置；

(i)判断是否编历模型节点，如果没有则继续执行步骤(h)，如果编历则继续下一步。