EQP：多项式时间步长函数的积分库及其应用

软件X 12（2020）100636原始软件出版物EQP-一个2D/3D库，用于多项式时间步长函数的积分Gregorio Mariggiès，Pastiano Fichera，Mauro CorradoPasta，Giulio Ventura意大利都灵理工大学结构、岩土和建筑工程系，地址：Corso Duca degli Abruzzi 24，10129 Torino，ar t i cl e i nf o文章历史记录：接收6十月2020收到修订版2020年11月14日接受2020年保留字：数值积分Heaviside阶跃函数多项式函数计算机图形库Fortran质量性质子域求积a b st ra ct本文提供的EQuivalent Polynomials库EQP是用于利用经典求积规则（例如，Gauss-Legendre）对由下式的乘积给出的函数进行数值积分的强大工具：任意多项式乘以海维赛德阶跃函数。该库可以处理多种形状在一维，二维和三维的积分域。最初由Ventura开发（Ventura，2006），以克服在扩展有限元法的背景下集成不连续函数的长期问题，EQP库最近已被推广以满足需求从计算力学到计算机图形学，几何区域（质量）属性的评估和一般的计算机模拟©2020作者由爱思唯尔公司出版这是CC BY-NC-ND下的开放获取文章许可证（http://creativecommons.org/licenses/by-nc-nd/4.0/）。代码元数据当前代码版本v1.2用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX-D-20-00060Code Ocean compute capsule无GNU通用公共许可证（GPL）使用的代码版本控制系统使用FORTRAN 90的编译要求、操作环境依赖性如果可用，链接到开发人员文档/手册www.equivalent-polynomials.net/reference-documentation问题支持电子邮件giulio. polito.it1. 动机和意义在许多物理和工程领域以及计算机图形学中，需要对多项式函数进行数值积分让我们考虑，作为几个代表性的例子，计算惯性矩和产品进入刚体的基于物理的动画所需的几何形状的惯性张量[1，2];计算力学系统的刚度矩阵的有限元方法的框架中的固体和结构的力学行为的预测[3，4];计算质量，总能量，角动量，和*通讯作者。电子邮件地址：mauro. polito.it（M. Corrado）。https://doi.org/10.1016/j.softx.2020.100636熵来施加控制大气动力学和热力学的守恒定律[5]。当定义域的形状是初等几何（三角形、直角三角形、平行六面体）或可以通过变换回到初等几何时，有效的数值求积规则可用于多项式被积函数。然而，当多项式显示跳跃不连续性或要在子域上积分时，引入乘法Heaviside阶跃函数。在这种情况下，由于数值求积规则（如例如，在计算机图形学中，当计算可以被看作规则几何形状的分区的复杂物体的几何属性时，可能会遇到这种情况让我们考虑一下，因为例如，图中的身体+。 1，它是通过分裂2352-7110/©2020作者。由爱思唯尔公司出版。这是一篇开放获取的文章，使用CC BY-NC-ND许可证（http://creativecommons.org/licenses/by-nc-nd/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softxG. Mariggiant，S.菲切拉湾Corrado等人软件X 12（2020）1006362Ω==+={∈H（x）=∈Q∑<$+pn（X）d=H（X）pn（X）d（7）Ω¯Ω公司简介Fig. 1. 由不连续表面Γ穿过的3D域。三次域由不连续面Γ分成两部分。可以通过积分来计算Σ+的几何性质，例如体积、力矩和惯性积：I=+pn（x）d=H（x）pn（x）d（1）其中pn（x）是关于计算所需几何性质（pn（x））所需的变量x（x， y， z）集合的n次多项式函数1表示体积，pn（x）y2z2对于相对于x轴的二次面积矩等），并且H（x）是标准Heaviside阶跃函数，由等式（1）给出（二）、由于H∞（x）pn（x）是在整个区域上连续的多项式函数，因此可以使用适当的数值求积规则对其进行精确积分[33]。注意，一般情况下，equiv-多项式函数H（x）与pn（x）具有相同的次数，因此在Eq. （6）与Eq. （1）它的数量增加了一倍。注意，引入等价多项式允许在标准域上积分，而不是在非标准分区子域上积分。这是等效多项式方法的优点。本文介绍了一个独立的库，EQP，它提供了，为五种不同形状的积分域，即三角形，圆形，四面体和paral-leleptape，等效多项式函数的表达式，H（x），作为不连续线位置的函数，或surface，r.关于确定等价多项式表达式所采用的数学程序的细节可以在参考文献中找到。[6、7]。2. 软件描述该库是用Fortran构建的，Fortran是一种多平台语言，可以与现有的库耦合。尽管如此，该库的结构是简单的，因此可以毫不费力地将其移植到任何语言。1if x+0其他智慧（二）2.1. 软件功能定义域的正部分由n+向量定义，它指向n+。矢量n+与内部不连续表面正交，并且在3D情况下，其分量由表面Γ的方程的系数a、b和c（3））。ax+by+cz+d=0（3）上述数值积分问题通常通过划分积分域以生成正交子单元来解决，其中被积函数是多项式。怎么--虽然EQP库的目的是提供等价多项式函数H∞（x）的表达式，但对于实际使用库需要被包装成一个算法，评估- uates在方程中定义的积分（6）通过应用例如高斯求积规则。尽管如此，该库可以直接包含在任何能够积分纯（不涉及不连续性）多项式函数的求积规则的算法EQP库的功能在本文中通过一般示例来呈现。假设多项式pn（X）为在子域上积分，子域通过分裂paral-然而，很难识别和分类所有可能的如图所示，有一个平面的lelepeppelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpelpel凌晨2X=积分子域的形状，特别是对于3D几何形状。因此，在[6，7]中，在以下的上下文中，已经提出了消除正交域的细分的方法：（X， Y， Z）是定义问题的全局参考系。所需的解决方案可以通过应用Eq.（6）、一旦已知等价多项式函数H（X）扩展有限元法（XFEM）[8，9]。该方法在其他框架中也被积极使用或讨论，一些非详尽的参考文献[10虽然I=Ω∫˜¯方法已开发的多项式积分，在原则上它可以应用/扩展到样条的分段多项式表示，用于引入不连续或修剪域。在这一领域，可以回顾关于构造特殊化求积规则的文献[28事实上，所提出的方法是独立的特定的数值积分中采用该方法是基于取代海维赛德函数H（x）与等价的多项式函数H（x），使得：首先，问题要转化为一个求积问题，一个标准域。因此，应用从全局参考系统XYZ到父坐标系统xyz的变量的改变，以便评估标准规则几何形状上的积分，如图所示。2b.事实上，这样的变换允许用单个亲本几何学来处理各种情况。例如，具有任何大小和位置的平行六面体可以被带回到在局部参考系（x， y， z）中定义的立方父几何，[-1，+1]如图所示。2 b.HΩH（x）pn（x）d联系我们 H（x）pn（x）d（4）库中实现变量变化的数学方法是常规的等参映射等价多项式H（x）依赖于Γ，具有相同的用有限元法求解[34]。 设P（x，y，z）为父坐标系中对应于pn（x）的次数，并由多项式函数表示：点P在全局坐标系中。仿射H（x）=c·m（x）（5）其中m（x）收集单项式基，例如m（x）=（1，x，y，z，x2，任何P（x，y，z）到P¯（X，Y，Z）的映射定义为：X=∑N（x，y，z）X（8a）. . . ），并且c是系数向量。考虑到Eq。（4），Eq.（1）可以改写为：我我i=q1I=H（x）pn（x）d（6）ΩY=Ni（x，y，z）Yi（8b）i=1Ω−G. Mariggiant，S.菲切拉湾Corrado等人软件X 12（2020）1006363∑∑||==||=必须正确定义，以便使单位向量n+指向期望的子域方向。例如，可以通过简单地更改所有不连续系数的符号此外，请注意，不连续性不一定要相交于<$$>：当不发生这种情况时，如果<$$>+=<$<$，则求积结果将为零或者，如果Σ+=Σ，则与Σ上的积分一致。图二. 平行六面体元素的映射：（a）全局坐标系中的配置;（b）父坐标系中的配置QZ=Ni（x，y，z）Zi（8c）i=1其中q表示在全局参考系中具有坐标（Xi，Yi，Zi）的几何元素的节点（顶点）的数量，并且Ni（x， y， z）表示根据父元素的第i个节点的局部坐标的插值函数类似地，不连续性方程Γ<$（X），定义为由用户在全局参考系中定义的，必须转换为在父坐标系中定义的Γ（x）（见图2）。在图书馆里，这是通过以下方式完成的：计算在所述全局参考系中的所述积分域的节点与所述积分域的节点之间的所述有符号距离Di用Di表示Γ<$（X）的系数a、b、c和d通过求解线性方程组;将Γ'（X）的变量X、Y和Z替换为ex-Eq. （8），从而获得由依赖于Di的新系数a′、b′、c′和d′定义的Γ（x）。例如，在三角形的情况下，可以得到：a′=D2−D1（9a）b′=D3−D1（9b）c′=D1（9c）在将Γ<$（X）变换为Γ（x）之后，由库提供等价多项式函数H<$（x）的适当表达式参照父域坐标系。然后，通过雅可比矩阵，在正交中引入坐标和积分域变换，雅可比矩阵包含插值函数Ni的偏导数，关于父系统变量微分x，y， z[35]：这是可以从参考文献中推断出来的。[6，7]中，可以与所提出的方法精确集成的多项式函数的次数和组成取决于为确定等效多项式而施加的一些条件表1中列出了包含在当前版本库中的父几何形状和可以组成要积分的多项式函数的单项式。然而，请注意，可以将库扩展到每个检查域中的任何多项式次数。2.2. 软件构架软件库的体系结构很简单。主库文件是eqpol.f90。它包含了将不连续性方程从全局坐标变换到父坐标的算法系统，以获得系数并计算等价多项式函数H∞（x）。它是由其他文件完成的，定义等效多项式函数所需的系数的代数表达式。为 了 演 示 使 用 ， 该 库 由 主 程 序 文 件 main.f90 和 模 块 文 件mapping_module.f90补充，它们不是该库的一部分，但允许计算方程中的积分（10）在某些情况下。图书馆的实际使用预计采取以下步骤：1. 初步数据准备：积分域的选择（见表1）;全局参考系中域节点坐标的个性化，或圆形域的中心和半径;在全局坐标参考中的不连续平面系数的个性化2. 变换到父（标准）域，并通过子例程Heqpol_coefficients评估等效多项式系数向量3. 正交通过任何选定的规则，例如，当量（10）中，由函数HeqPol提供在正交点处的等价多项式的值，由函数det_J给出雅可比行列式的变换行列式的值。该库的当前版本允许在以下假设下评估精确的求积结果：IH<$（X）pn（X）d<$H<$（x）pn（x）J d<$Ω Ω中文（简体）转换雅可比矩阵在方程中。（10）是常数（这总是适用于域：三角形、圆形、正方形、四面体、平行六面体）;wjH<$（xj，yj，zj）pn（xj，yj，zj）|J（xj，yj，zj）|j=1其中J是雅可比行列式的行列式。方程的积分（7）在Eq. （10）用标准形式的由方程式其中，gp表示Gauss-Legendre正交点的数量值得再次强调的是，该库在整个域上积分，并产生子域上积分的结果。因此，不连续性方程多项式pn是表1中列出的单项式的线性组合。注意，如果用户想要在他自己的正交算法中直接集成库的核心 ， 则 必 须 预 见 两 个 调 用 ： 对 每 个 积 分 域 的 子 例 程Heqpol_coefficients的调用（上面的步骤2），以及对每个正交点的函数HeqPol的另一个调用（上面的步骤3）。函数det_J不包含在库的核心中，因为它应该是正交算法的一部分。········G. Mariggiant，S.菲切拉湾Corrado等人软件X 12（2020）1006364<$=表1积分整环，整环类型，父整环，单项基。图三. 坐标图（etype 21）：（a）全局坐标系中的配置;（b）父坐标系中的配置3. 说明性实例以下示例基于库用法示例文件main.f90。对于所提供的示例，库提供了机器精度的精确结果。3.1. 平行四边在该示例中，图1A中所示的图1B中所示的图1C中所示的图1B中所示的3（a）是一个考虑。 尺寸为米。 元素被方程4的直线X+Y−26=0变成了Y+和Y−。 目标是计算面积和惯性张量。平行四边形被软件映射到父坐标系上，如图所示。 3（b），其中H（x）被计算，积分进行。启动库示例程序后，通过在屏幕上键入以下内容提供输入等价多项式库etype（20，21，22，30，31）：21a，b，c：4，1，-26按照如下所示的方案插入元素坐标：4----------3||||||1----------2X（1），Y（1）：4，4X（2），Y（2）：7，5X（3），Y（3）：8，8X（4），Y（4）：5，7在此阶段，用户可以选择对表1中列出的所有单项式进行积分， 用于所选 几何类型， 将出现 在屏幕上的 单项式， 或在user_fun.txt文件中定义为表1中列出的单项式的线性组合的多项式。例如，如果用户想要计算区域的面积A，则单项式p应输入1在user_fun.txt文件中，而单项式y2，x2和xy具有来得到惯性张量I。相对于XY轴计算的面积和惯性张量并由图书馆输出报告的是：其中，数字21标识了数字，而数字21标识了数字。A= 6。418米2（11）Bers4，1，-26是直线的系数，I=ρ[254. 436-254。281]（单位：千克平方米）（12）全球坐标系之后，库要求输入全局坐标系中积分域的节点坐标：S -254281258. 876其中ρs表示材料的表面密度（kg/ m2）G. Mariggiant，S.菲切拉湾Corrado等人软件X 12（2020）1006365- − − +=[客户端]图四、 被一个平面分裂的立方体：（a）正部分，+，积分域;（b）负部分，−，不被积分。3.2. 平行六面体本文分析了一个3D示例：边2 m的立方体被切割为方程XYZ5的曲面。5、如图所示。 四、输入过程类似于前面的示例，如下所示等价多项式库etype（20，21，22，30，31）：31a、b、c、d：-1、-1、-1、+5.5按照如下所示的方案插入元素坐标：5-----------8/|/|//下一页|//下一页|6-----------7||1--------|--4| //下一页|//下一页|/|/2-----------3X（1）、Y（1）、Z（1）：3,1,0X（2）、Y（2）、Z（2）：3,3,0X（3）、Y（3）、Z（3）：1,3,0X（4）、Y（4）、Z（4）：1,1,0X（5）、Y（5）、Z（5）：3,1,2X（6）、Y（6）、Z（6）：3,3,2X（7）、Y（7）、Z（7）：1,3,2X（8）、Y（8）、Z（8）：1,1,2的体积V和惯性张量I（图1）。 4（a））相对于XYZ轴计算并由库输出报告的是：V= 5。458米3（13）这里设想的应用涉及计算几何学的背景，其中EQP库可用于计算复杂图形的几何特性，所述复杂图形通过用平面切割基本形状（例如正方形或立方体）并最终将它们组合在一起而获得。此外，库的特性可以在涉及对象的形状和位置的动态变化的模拟中被利用，例如，将对象或离散化单元分解成几块。在这种情况下，例如，已经有使用XFEM在3D元素中模拟脆性断裂和表面裂纹模式而无需重新网格划分的经验[37]，优化了[38]中提出的方法。因此，EQP库可以集成在许多现代计算工具，广泛的应用领域。5. 结论本文提出的EQP库是一个现成的工具，大大简化了通过分裂一个标准的正交域得到的子域上的多项式函数的积分的数值计算，通过将它们恢复到整个域上的等价多项式的积分。EQP的优点在于计算速度快、不需要复杂的子域计算、数学方法的绝对通用性和可扩展性，有望成为一种广泛应用的通用工具。该软件的未来发展设想通过添加新的改进算法来自动过滤，校准，映射和组合现有的几何形状，以轻松创建多体复杂的几何形状。竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作致谢感谢意大利教育、大学和研究部提供的研究资助“XFAST-SIMS：复杂结构系统的快速准确模拟”。感谢意大利都灵理工大学通过分散研究补助金倡议提供的研究支持引用[1]哈恩JK。刚体的真实动画。比较图表1988;22：299 - 308。[2]Bender J，Erleben K，Trinkle J.计算机图形学中刚体动力学的交互式模拟。Comput Graph Forum2014;33：246-70.I=ρ二十五008-17。791-7。838-17。79125. 008-7。838-7。838- 7。83839. 456G. Mariggiant，S.菲切拉湾Corrado等人软件X 12（2020）1006366（单位：千克平方米）（14）[3]给K-J洗澡有限元程序。第二版Watertown，MA; 2014.[4]Belytschko T，Liu WK，Moran B，Elkhovik K.连续体和结构的非线性有限元。第二版Wiley; 2014.[5]张文，张文，等.大气物理学式中，ρ表示材料的体积密度（kg/m3）。4. 影响EQP库是一个有用的工具，积分不连续函数与任何数值求积方法 ， 而 不 分 裂 的 积 分 域 。 它 的 效 率 已 经 在 计 算 力 学（XFEM/GFEM）领域得到了证明[6，7]。然而，它的影响要广泛得多，因为多项式函数的数值积分是许多领域中的常见问题的rology 。 In： Temam RM ， Tribbia JJ， editors. 大 气 和 海 洋 的 计 算 方 法 。Elsevier Science; 2008，p.1-120[6] 文图拉湾关于不连续的正交子胞的消除，扩展有限元法中的一些函数。Internat J NumerMethods Engrg2006;66：761-95.[7]范图拉G，本韦努蒂E. Heaviside函数求积的等价多项式. Internat J NumerMethods Engrg2015;102：688-710.[8]Belytschko T，Black T.最小网格重划有限元中的弹性裂纹扩展。国际数值方法工程杂志1999;45（5）：601-20.[9]张文辉，王文辉，王文辉.材料建模的扩展/广义有限元方法综述。模型模拟材料科学工程2009;17：1-24.G. Mariggiant，S.菲切拉湾Corrado等人软件X 12（2020）1006367[10]吕建辉，焦玉英，拉布丘克，庄晓英，冯晓堂，谭芳。PU基数值方法中三维裂纹奇异性数值积分的通用算法。第363章：你是我的女人112908-1[11]Düster A，Allix O.力矩拟合的选择性丰富及其应用来切割有限元和单元。Comput Mech2020;65（2）：429-50.[12][10] Ali T，Mostefa B，Abdelkader D，Abdelkrim A，Habibe K.使用XFEM对铝板进行实验和数值断裂建模。Int J Eng ResAfrica2020;46：45-52.[13]Surendran M，Natarajan S，Palani G，Bordas S.模拟疲劳裂纹扩展的线性光滑扩展有限元法。 工程机械2019;206：551-64。[14]Müller B，Krämer-Eis S，Kummer F，Oberlack M.可压缩流体浸没边界的高阶间断Galerkin方法。Internat J Numer Methods Engrg2017;110（1）：3-30.[15]萨 耶 河 超 矩 形 中 隐 式 定 义 曲 面 和 体 积 的 高 阶 求 积 方 法 。 SIAM J SciComput2015;37（2）：A993-1019.[16]张文，张文，张文，等.固体力学中X-有限元积分格式的应用.北京：机械工业出版社，2001。计算方法应用机械工程2015;283：551-72.[17]吴晓波，李晓波，李晓波.一种精确、鲁棒、易于实现的任意多面体积分方法：应用于嵌入式接口方法。J Comput Phys2014;273：393-415.[18]本 韦 努 蒂 湖 基 体 - 夹 杂 物 界 面 等 效 本 征 应 变 扩 展 有 限 元 法 。 ComputMech2014;53（5）：893-908.[19]Alves P，Barros F，Pitangueira R.广义有限元法的面向对象方法。Adv EngSoftw2013;59：1-18.[20]Loehnert S，Mueller-Hoeppe D，Wriggers P.三维修正XFEM方法和有限变形理论的扩展。Internat J Numer MethodsEngrg2011;86（4-5）：431-52.[21]Mousavi S，Sukumar N.不规则凸多边形和多面体上多项式和间断函数的数值积分。ComputMech2011;47（5）：535-54.[22][10] Antoietti P，Verani M，Vergara C，Zonca S.用多边形网格上的高阶间断Galerkin方法数值求解流固耦合问题。有限的元素肛门Des 2019;159：1[23][10]李文辉，李文辉.线弹性断裂力学的离散和相场方法：比较研究和最新评论。应用科学（瑞士）2019;9（12）。[24]吴建勇，邱建芳，阮伟，曼达尔，庄立杰。固体局部破坏的数值模拟：XFEM与PF-CZM。 计算方法应用机 械 工 程 2019;345：618-43。[25]Formaggia L，Vergara C，Zonca S.复合材料网格的非拟合扩展有限元。Comput Math Appl2018;76（4）：893-904.[26]张伟，王伟，王伟.扩展有限元法中无单元划分的裂纹不连续性建模。InternatJNumer Methods Engrg2017;110（11）：1021-48.[27]Künner L，Zander N，Kollmannsberger S，Rank E.智能八叉树：在3D中精确集成不连续函数。计算方法应用机 械 工 程 2016;306：406-26.[28]Hiemstra RR，Calabro F，Schillinger D，Hughes TJ. 等几 何 分 析 中 张量积和分级细化样条的最佳和简化求积规则。计算方法应用机械工程2017;316：966-1004。[29]Bartoregon M，Calo VM.二次和三次样条空间的Gauss-Galerkin求积规则及其在等几何分析中的应用。2017年12月22日，《明史》卷57。[30]Bartoregon M，Calo VM.奇次样条空间的最佳求积规则及其在基于张量积的等几何分析中的应用。 计算方法应用机 械 工 程 2016;305：217-40.[31]放大图片作者：Calabro F，Sangalli G，Tani M.加权求积法快速构造等几何伽辽金矩阵。计算方法应用机械工程2017;316：606-22.[32]Scholz F，Mantzaflaris A，Jüttler B. 等几何分析中截断求积的一阶误差修正。在：开姆尼茨精细元素研讨会. Springer; 2017，p. 297-321[33]斯特劳德河重积分的近似计算。Englewood Cliffs，NJ：Prentice Hall;1971.[34]Zienkiewicz OC，Taylor RL，Zhu JZ.有限元法：其基础和基本原理。第六版Butterwest-Heinemann;2005.[35]作者：J.高等微积分及其在工程和物理科学中的应用。John Wiley Sons Inc;1980.[36]Davis PJ，Rabinowitz P.数值积分方法。第二版中国科学院出版社，1984年.[37]宋C，张H，吴Y，包H.无需重新网格化的切割和压裂模型。In：Chen F，Jüttler B，editors.几何建模与处理进展。Springer Berlin Heidelberg; 2008，p.107比18[38]Iben HN，O'Brien JF.产生表面裂纹图案。Graph Models 2009;71（6）：198-208. http://dx.doi.org/10.1016/j.gmod.2008.12.005网站。