CN106202689A - 一种软组织有限元模型的加速计算方法 - Google Patents
一种软组织有限元模型的加速计算方法 Download PDFInfo
- Publication number
- CN106202689A CN106202689A CN201610524437.4A CN201610524437A CN106202689A CN 106202689 A CN106202689 A CN 106202689A CN 201610524437 A CN201610524437 A CN 201610524437A CN 106202689 A CN106202689 A CN 106202689A
- Authority
- CN
- China
- Prior art keywords
- soft tissue
- sequence number
- unit
- model
- array
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/06—Power analysis or power optimisation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明涉及一种人体组织有限元仿真技术，特别涉及一种软组织有限元模型的加速计算方法领域。本发明包括以下步骤：建立软组织网格模型的四面体单元和三角形单元的对应关系；进行参与计算单元的筛选；建立参与计算的单元的新四面体单元序号、新节点序号与在原模型中对应的四面体单元序号、节点序号的关系；建立平衡方程，并进行形变计算。本发明能够有效地减少软组织模型的计算规模和计算时间，同时也能够真实地仿真出人体软组织的生物力学特性，对虚拟手术培训系统的建立具有一定的意义。
Description
技术领域
本发明涉及一种人体组织有限元仿真技术，特别涉及一种软组织有限元模型的加速计算方法领域。
背景技术
虚拟手术培训系统能够提高医生的手术技能，从而降低手术失败的风险。构建一个可靠的虚拟手术培训系统，不仅需要精确地反映人体组织真实的力学特性，还有保证其仿真的实时性。人体软组织具有复杂的力学特性，而利用有限元法能够精确地预测人体组织的变形，但由于其较长的计算时间，难以在虚拟手术培训系统中直接应用。
对现有技术的文献检索发现，闫桢南在上海交通大学的硕士学位论文《软组织非线性形变快速模型与应用》中，基于有限元法将软组织模型分为操作区域和非操作区域，并只对操作区域的形变进行计算的方法来减少软组织模型计算规模，从而达到加速仿真的效果。而这种方法需要用户在模拟手术演练前，根据自己以往的经验，对软组织的模型进行手动划分操作区域，因此缺乏灵活性，尤其是对手术路径较长的情况，其加速效果并不理想。
综上所述，目前在建立虚拟手术培训系统的过程中，软组织的实时仿真仍然是一项巨大的挑战。
发明内容
本发明的目的在于为了克服现有技术中的不足，提供一种软组织有限元模型的加速计算方法，能够同时实现虚拟手术培训系统的真实性和实时性。
本发明为解决上述技术问题采取的技术方案是：
一种软组织有限元模型的加速计算方法，其特征在于包括以下步骤：
步骤1）：根据软组织的网格模型，创建数组A，将软组织模型的每个四面体单元序号与每个四面体单元所包含的四个三角形单元序号相连，使得任意一个四面体单元序号可得出其包含的四个三角形单元序号。同时创建数组B，将软组织模型的每个三角形单元序号与共用每个三角形单元的两个相邻四面体单元的序号相连，使得通过任意一个三角形单元序号可得出共用其三角形单元的两个相邻四面体单元序号，如果其中某一三角形单元为软组织模型的边界单元，则只与包含当前三角形单元的四面体单元序号相连。
步骤2）：筛选参与计算的单元的集合，以受力点所在单元的内心为基准点，给出如下所示的筛选条件的关系式：
(1)
式中，l为当前四面体单元的内心与受力点所在单元的内心之间的距离，a为筛选系数，R为受力单元的外接圆半径；
判断当前的四面体单元是否满足筛选条件关系式（1），如果是，将当前四面体单元列入参与计算单元的集合，判断首先从受力点所在四面体单元的相邻四面体单元开始进行，然后再对每个相邻四面体单元相邻的四面体单元进行判断，参与计算单元集合的判断以受力点所在四面体单元为中心，由里到外进行，直到不满足筛选条件关系式（1）为止。筛选参与计算的单元时，根据筛选判断的先后顺序，将参与计算的单元集合的四面体单元序号和节点序号进行重新排序。当受力点所在四面体单元个数大于1时，依次进行受力点所在四面体单元个数次的循环进行判断，并统一封装到参与计算单元的集合；每个四面体单元的相邻四面体单元，通过步骤1）中生成的数组A和数组B进行查找。
步骤3）：创建数组C，将步骤2）中生成的参与计算单元集合的新四面体单元序号和原模型中的四面体单元序号相连。创建数组D，将步骤2）中生成的参与计算单元集合的新节点序号和原模型中的节点序号相连。
步骤4）：基于步骤2）所生成的参与计算的单元集合，并结合步骤3）中所生成的数组C和数组D，封装整体刚度矩阵，在施加的外加载荷和边界条件下，只对软组织模型中所筛选出的参与计算单元的集合进行形变计算。
上述的一种软组织有限元模型的加速计算方法，所述步骤1）中的数组A为m行4列的二维数组，数据类型为int型，m为软组织模型所包含的四面网格单元的个数，每个行下标代表软组织模型的四面体网格单元的序号，每行的四个列下标所对应的元素分别存储每个行下标所对应四面体单元所包含的四个三角形单元的序号。数组B为n行2列的二维数组，数据类型为int型，n为软组织模型所包含的三角形单元的个数，每个行下标代表软组织模型的三角形网格单元的序号，每行的两个列下标所对应的元素分别存储共用当前三角形单元的两个四面体单元的序号，如果当前三角形单元为软组织模型边界面上的单元，则第一个元素存储包含当前三角形单元的四面体单元的序号，并且第二个元素为-1。
上述的一种软组织有限元模型的加速计算方法，所述步骤3)中的数组C为int型的一维数组，数组C的下标表示原模型的四面体单元序号，每个参与计算单元的新四面体单元序号分别存储到对应的原模型的四面体单元序号下标的数组C的元素中，而原模型的四面体单元序号中，非参与计算单元下标的元素值赋为-1。数组D为int型的一维数组，数组D的下标表示原模型的节点序号，每个参与计算单元的新节点序号分别存储到对应的原模型节点序号下标的数组D的元素中，而原模型节点序号中，非参与计算节点下标的对应元素值赋为-1。
上述的一种软组织有限元模型的加速计算方法，所述步骤4）中软组织形变计算的过程中，软组织的模型采用指数多项式混合形式的超弹性模型，模型的表达式如下：
(2)
式中，W为应变能密度函数，C1、β为模型的材料常数，I1和I2为应变张量的主不变量；
基于软组织的指数多项式混合形式的超弹性模型建立平衡方程，然后施加边界条件和外加载荷，可进行形变计算。在进行形变计算的过程中，平衡方程采用完全拉格朗日法，进行分步载荷求解，每一载荷步的非线性方程组采用修正的牛顿迭代法进行求解。在软组织有限元计算的过程中，整体刚度矩阵采用GPU计算，非线性方程组采用CPU进行求解。
本发明的有益效果是：
本发明中，基于指数多项式混合形式的超弹性模型，进行人体软组织的有限元仿真，能够在一定程度上保证其仿真精度。然后，根据圣维南原理，通过给出一个筛选系数a来筛选软组织模型中参与计算单元的集合，并只对局部筛选出的参与计算的单元集合进行有限元计算，
提高软组织有限元模型的计算速度。用户通过合理地选取筛选系数a的值，能够有效地减少系统的计算规模，同时能够真实地仿真出人体软组织的生物力学特性。整个过程中，整体刚度矩阵采用GPU进行计算，能够进一步地提高了整个有限元计算的速度。
附图说明
图1是本发明流程框图，即软组织加速计算流程图；
图2是筛选参与计算单元集合的判断示意图；
图3是本发明应用的平面示意图。
具体实施方式
具体实施方式一：如图1所示，本实施方式所述的一种软组织有限元模型的加速计算方法，包括以下步骤：
步骤1）：根据软组织的网格模型中的每个四面体单元序号与三角形单元序号，创建数组A，将软组织模型的每个四面体单元序号与每个四面体单元所包含的四个三角形单元序号相连，使得通过软组织模型的任意一个四面体单元序号可得出其包含的四个三角形单元序号。同时创建数组B，将软组织模型的每个三角形单元序号与共用每个三角形单元的两个相邻四面体单元序号相连，使得通过软组织模型的任意一个三角形单元序号可得出共用该三角形单元的两个相邻四面体单元序号，如果其中某一三角形单元为软组织模型的边界单元，则只与包含当前三角形单元的四面体单元序号相连。本步骤在离线状态下实现，因此对软组织计算速度无影响。
步骤2）：筛选参与计算的单元的集合，以受力点所在单元的内心为基准点，给出如下所示的筛选条件的关系式：
(1)
式中，l为当前四面体单元的内心与受力点所在单元的内心之间的距离，a为筛选系数，R为受力单元的外接圆半径。
判断当前的四面体单元是否满足筛选条件关系式（1），如果是，将当前四面体单元列入参与计算单元的集合，判断首先从受力点所在四面体单元的相邻四面体单元开始进行，然后再对每个相邻四面体单元相邻的四面体单元进行判断。如图2所示，受力点所在的四面体单元包含i、j、k和p四个点，并与相邻的四面体单元共用由i、k和p组成的三角形单元，点n1到n2的距离即为公式(1)中的，并且满足时，可判断该相邻四面体单元为参与计算的单元。参与计算单元集合的判断以受力点所在四面体单元为中心，由里到外进行，直到不再满足筛选条件关系式（1）为止。筛选参与计算的单元时，根据筛选判断的先后顺序，将参与计算的单元集合的四面体单元序号和节点序号进行重新排序。当受力点所在四面体单元个数大于1时，依次进行受力点所在四面体单元个数次的循环进行判断，并统一封装到参与计算单元的集合。每个四面体单元的相邻四面体单元，通过步骤1）中生成的数组A和数组B进行查找，即通过数组A寻找当前四面单元所包含的四个三角形面单元，再根据得出的四个三角形面单元寻找与当前面四面体单元相邻的四个四面体单元，从而实现相邻单元的快速查找，并依次进行筛选条件的判断。
步骤3）：创建数组C，将步骤2）中生成的参与计算单元集合的新四面体单元序号和原模型中的四面体单元序号相连，创建数组D，将步骤2）中生成的参与计算单元集合的新节点序号和原模型中的节点序号相连，能够实现快速查询筛选出的参与计算的四面体单元新序号和节点新序号与原模型对应的原序号，并能够保证实时更新发生形变的节点坐标。
步骤4）：基于步骤2）所生成的参与计算的单元集合，并结合步骤3）中所生成的数组C和数组D，封装整体刚度矩阵，实现根据施加的外加载荷和边界条件，只对软组织模型中所筛选出的参与计算单元的集合进行形变计算。
具体实施方式二：本实施方式所述的一种软组织有限元模型的加速计算方法，所述的步骤1）中的数组A为m行4列的二维数组，数据类型为int型，m为软组织模型所包含的四面网格单元的个数，每个行下标代表软组织模型的四面体网格单元的序号，每行的四个列下标所对应的元素分别存储行下标对应四面体单元所包含的四个三角形单元的序号。例如，A[2][2]存储序号为2（假设四面体单元序号从0开始）的四面体单元中第三个三角形单元的序号。数组B为n行2列的二维数组，数据类型为int型，n为软组织模型所包含的三角形单元的个数，每个行下标代表软组织模型的三角形网格单元的序号，每行的两个列下标所对应的元素分别存储共用当前三角形单元的两个四面体单元的序号，如果当前三角形单元为软组织模型边界面上的单元，则第一个元素存储包含当前三角形单元的四面体单元的序号，并且第二个元素为0。例如，B[2][1]存储共用序号为2（假设三角形单元序号从0开始）的三角形单元的两个相邻四面体单元中，第2个四面体单元的序号，如果该三角形单元为软组织模型的边界单元，则B[2][1]存储的值为-1。
具体实施方式三：本实施方式所述的一种软组织有限元模型的加速计算方法，所述步骤3)中的数组C为int型的一维数组，数组C的下标表示原模型的四面体单元序号，每个参与计算单元的新四面体单元序号分别存储到对应的原模型四面体单元序号下标的数组C的元素中，而原模型四面体单元序号中，非参与计算的单元下标的对应元素值赋为-1。例如，参与计算单元集合中序号为2的四面体单元在原模型中的序号为7（假设四面体单元序号从0开始），则C[7]赋值为2，如果原模型中序号为7的四面体单元不是参与计算单元集合中的四面体单元，则C[7]赋值为0。数组D为int型的一维数组，数组D的下标表示原模型的节点序号，每个参与计算单元的新节点序号分别存储到对应的原模型节点序号下标的数组D的元素中，而原模型节点序号中，非参与计算节点下标的元素值赋为0。例如，参与计算单元集合中序号为2的节点在原模型中的序号为7（假设节点序号从0开始），则D[7]赋值为2，如果原模型中序号为7的节点不在参与计算单元集合中的四面体单元上，则D[7]赋值为-1。
具体实施方式四：本实施方式所述的一种软组织有限元模型的加速计算方法，所述步骤4）中，如图3所示，只对软组织的有限元模型中参与计算单元集合，即图3中的灰色区域进行计算。软组织的模型采用指数多项式混合形式的超弹性模型，模型的表达式如下：
(2)
式中，W为应变能密度函数，C1、β为模型的材料常数，I1和I2为应变张量的主不变量；
基于软组织的指数多项式混合形式的超弹性模型建立平衡方程，然后施加边界条件和外加载荷，可进行形变计算。在进行形变计算的过程中，平衡方程采用完全拉格朗日法，进行分步载荷求解，每一载荷步的非线性方程组采用修正的牛顿迭代法进行求解，而在求解非线性方程组的过程中产生的线性方程组可采用CG法进行求解。在软组织有限元计算的过程中，整体刚度矩阵采用GPU计算，非线性方程组采用CPU进行求解，来进行进一步地提高计算效率。
Claims (4)
1.一种软组织有限元模型的加速计算方法，其特征在于包括以下步骤：
步骤1）：根据软组织的网格模型，创建数组A，将软组织模型的每个四面体单元序号与每个四面体单元所包含的四个三角形单元序号相连，使得任意一个四面体单元序号可得出其包含的四个三角形单元序号；同时创建数组B，将软组织模型的每个三角形单元序号与共用每个三角形单元的两个相邻四面体单元的序号相连，使得通过任意一个三角形单元序号可得出共用其三角形单元的两个相邻四面体单元序号，如果其中某一三角形单元为软组织模型的边界单元，则只与包含当前三角形单元的四面体单元序号相连；
步骤2）：筛选参与计算的单元的集合，以受力点所在单元的内心为基准点，给出如下所示的筛选条件的关系式：
(1)
式中，l为当前四面体单元的内心与受力点所在单元的内心之间的距离，a为筛选系数，R为受力单元的外接圆半径；
判断当前的四面体单元是否满足筛选条件关系式（1），如果是，将当前四面体单元列入参与计算单元的集合，判断首先从受力点所在四面体单元的相邻四面体单元开始进行，然后再对每个相邻四面体单元相邻的四面体单元进行判断，参与计算单元集合的判断以受力点所在四面体单元为中心，由里到外进行，直到不满足筛选条件关系式（1）为止；筛选参与计算的单元时，根据筛选判断的先后顺序，将参与计算的单元集合的四面体单元序号和节点序号进行重新排序；当受力点所在四面体单元个数大于1时，依次进行受力点所在四面体单元个数次的循环进行判断，并统一封装到参与计算单元的集合；每个四面体单元的相邻四面体单元，通过步骤1）中生成的数组A和数组B进行查找；
步骤3）：创建数组C，将步骤2）中生成的参与计算单元集合的新四面体单元序号和原模型中的四面体单元序号相连；创建数组D，将步骤2）中生成的参与计算单元集合的新节点序号和原模型中的节点序号相连；
步骤4）：基于步骤2）所生成的参与计算的单元集合，并结合步骤3）中所生成的数组C和数组D，封装整体刚度矩阵，在施加的外加载荷和边界条件下，只对软组织模型中所筛选出的参与计算单元的集合进行形变计算。
2.如权利要求书1所述的一种软组织有限元模型的加速计算方法，其特征在于：所述步骤1）中的数组A为m行4列的二维数组，数据类型为int型，m为软组织模型所包含的四面网格单元的个数，每个行下标代表软组织模型的四面体网格单元的序号，每行的四个列下标所对应的元素分别存储每个行下标所对应四面体单元所包含的四个三角形单元的序号；数组B为n行2列的二维数组，数据类型为int型，n为软组织模型所包含的三角形单元的个数，每个行下标代表软组织模型的三角形网格单元的序号，每行的两个列下标所对应的元素分别存储共用当前三角形单元的两个四面体单元的序号，如果当前三角形单元为软组织模型边界面上的单元，则第一个元素存储包含当前三角形单元的四面体单元的序号，并且第二个元素为-1。
3.如权利要求书2所述的一种软组织有限元模型的加速计算方法，其特征在于：所述步骤3)中的数组C为int型的一维数组，数组C的下标表示原模型的四面体单元序号，每个参与计算单元的新四面体单元序号分别存储到对应的原模型的四面体单元序号下标的数组C的元素中，而原模型的四面体单元序号中，非参与计算单元下标的元素值赋为-1；数组D为int型的一维数组，数组D的下标表示原模型的节点序号，每个参与计算单元的新节点序号分别存储到对应的原模型节点序号下标的数组D的元素中，而原模型节点序号中，非参与计算节点下标的对应元素值赋为-1。
4.如权利要求书3所述的一种软组织有限元模型的加速计算方法，其特征在于：所述步骤4）中软组织形变计算的过程中，软组织的模型采用指数多项式混合形式的超弹性模型，模型的表达式如下：
(2)
式中，W为应变能密度函数，C1、β为模型的材料常数，I1和I2为应变张量的主不变量；
基于软组织的指数多项式混合形式的超弹性模型建立平衡方程，然后施加边界条件和外加载荷，可进行形变计算；在进行形变计算的过程中，平衡方程采用完全拉格朗日法，进行分步载荷求解，每一载荷步的非线性方程组采用修正的牛顿迭代法进行求解；在软组织有限元计算的过程中，整体刚度矩阵采用GPU计算，非线性方程组采用CPU进行求解。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610524437.4A CN106202689A (zh) | 2016-07-06 | 2016-07-06 | 一种软组织有限元模型的加速计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610524437.4A CN106202689A (zh) | 2016-07-06 | 2016-07-06 | 一种软组织有限元模型的加速计算方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN106202689A true CN106202689A (zh) | 2016-12-07 |
Family
ID=57466236
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610524437.4A Pending CN106202689A (zh) | 2016-07-06 | 2016-07-06 | 一种软组织有限元模型的加速计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106202689A (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109033742A (zh) * | 2018-06-21 | 2018-12-18 | 哈尔滨理工大学 | 一种用于模拟软组织剪切变形的超弹性模型 |
| CN110867218A (zh) * | 2018-08-27 | 2020-03-06 | 中国石油化工股份有限公司 | 一种分子电子能量计算方法及系统 |
| CN111090906A (zh) * | 2019-10-10 | 2020-05-01 | 惠州市德赛西威汽车电子股份有限公司 | 一种以非线性弹簧代替胶水实体建模的方法 |
-
2016
- 2016-07-06 CN CN201610524437.4A patent/CN106202689A/zh active Pending
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109033742A (zh) * | 2018-06-21 | 2018-12-18 | 哈尔滨理工大学 | 一种用于模拟软组织剪切变形的超弹性模型 |
| CN109033742B (zh) * | 2018-06-21 | 2019-06-21 | 哈尔滨理工大学 | 一种用于模拟软组织剪切变形的超弹性模型的建立方法 |
| CN110867218A (zh) * | 2018-08-27 | 2020-03-06 | 中国石油化工股份有限公司 | 一种分子电子能量计算方法及系统 |
| CN111090906A (zh) * | 2019-10-10 | 2020-05-01 | 惠州市德赛西威汽车电子股份有限公司 | 一种以非线性弹簧代替胶水实体建模的方法 |
| CN111090906B (zh) * | 2019-10-10 | 2023-07-07 | 惠州市德赛西威汽车电子股份有限公司 | 一种以非线性弹簧代替胶水实体建模的方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US7634394B2 (en) | Method of analysis of comfort for virtual prototyping system | |
| CN105302974B (zh) | 一种基于有限元和时变模态分析的柔性物体实时切割仿真方法 | |
| CN107590853A (zh) | 一种城市建筑群震害高真实度展示方法 | |
| CN106202689A (zh) | 一种软组织有限元模型的加速计算方法 | |
| CN105264533A (zh) | 用于心脏机电学的交互计算的方法和系统 | |
| CN102262699A (zh) | 基于无网格伽辽金与质点弹簧耦合的软组织形变仿真方法 | |
| Dzwinel et al. | Supermodeling in simulation of melanoma progression | |
| CN108511076B (zh) | 一种基于力学刺激和生物联合刺激的骨折愈合仿真系统 | |
| CN106777584B (zh) | 一种模拟骨折愈合过程的仿真系统 | |
| WO2005088485A2 (en) | Virtual prototyping system and method | |
| CN110289104B (zh) | 软组织按压和形变恢复的模拟方法 | |
| CN107133397B (zh) | 一种基于ale法对生物瓣膜进行双向流固耦合分析的方法 | |
| CN103699776A (zh) | 一种面向心血管介入手术仿真的导丝模拟方法 | |
| CN106777583B (zh) | 一种基于细胞活动的骨折愈合仿真方法 | |
| CN107134010A (zh) | 一种基于有限元的弹性软组织的形貌效果预测方法 | |
| CN104536831A (zh) | 一种基于多目标优化的多核SoC软件映射方法 | |
| CN101425188A (zh) | 带刚体核的广义弹簧振子形变仿真方法 | |
| CN114970381A (zh) | 运用cfd技术模拟左室附壁血栓的血流动力学特征的方法 | |
| CN103678888B (zh) | 一种基于欧拉流体模拟算法的心脏血液流动示意显示方法 | |
| DeCarlo et al. | Integrating Anatomy and Physiology for Behavior Modelling | |
| CN104463934A (zh) | 一种“质点-弹簧”系统驱动的点集模型动画自动生成方法 | |
| CN109829232A (zh) | 基于随机森林算法的分层布料模拟方法 | |
| CN105389853A (zh) | 一种基于多gpu的人脑变形仿真方法 | |
| CN106910233B (zh) | 一种虚拟昆虫动画角色的运动仿真方法 | |
| CN110047145A (zh) | 基于深度学习和有限元建模的组织变形模拟系统及方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20161207 |
|
| WD01 | Invention patent application deemed withdrawn after publication |