From 7657d65faf904eb963bd5f2d890258d43ce87f31 Mon Sep 17 00:00:00 2001 From: "Saullo G. P. Castro" Date: Wed, 11 Oct 2023 11:14:54 +0200 Subject: [PATCH 1/8] MAINT: newest Extension, classifiers --- setup.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/setup.py b/setup.py index 904db89..8ba222e 100644 --- a/setup.py +++ b/setup.py @@ -3,7 +3,7 @@ import inspect import subprocess from setuptools import setup, find_packages -from distutils.extension import Extension +from setuptools.extension import Extension from Cython.Build import cythonize @@ -80,6 +80,7 @@ def read(fname): Intended Audience :: Developers Intended Audience :: End Users/Desktop Topic :: Scientific/Engineering +Topic :: Scientific/Engineering :: Mathematics Topic :: Education Topic :: Software Development Topic :: Software Development :: Libraries :: Python Modules @@ -183,17 +184,19 @@ def read(fname): 'pyfe3d': ['*.pxd', '*.pyx'], '': ['tests/*.*'], } + keywords = [ - "finite elements", - "structural analysis", - "structural optimization", - "static analysis", - "buckling", - "vibration", - "structural dynamics", - "implicit time integration", - "explicit time integration", + 'finite elements', + 'structural analysis', + 'structural optimization', + 'static analysis', + 'buckling', + 'vibration', + 'structural dynamics', + 'implicit time integration', + 'explicit time integration', ] + s = setup( name = "pyfe3d", version = fullversion, From 92ece28820e960e211c103bb0e8ca684aff8cb06 Mon Sep 17 00:00:00 2001 From: "Saullo G. P. Castro" Date: Wed, 11 Oct 2023 23:16:17 +0200 Subject: [PATCH 2/8] DEV: piston theory for Quad4R element --- pyfe3d/quad4r.pyx | 2441 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 2436 insertions(+), 5 deletions(-) diff --git a/pyfe3d/quad4r.pyx b/pyfe3d/quad4r.pyx index bf4c4d1..f6df0e4 100644 --- a/pyfe3d/quad4r.pyx +++ b/pyfe3d/quad4r.pyx @@ -37,14 +37,30 @@ cdef class Quad4RData: M_SPARSE_SIZE, : int ``M_SPARSE_SIZE = 480`` + KA_BETA_SPARSE_SIZE, : int + ``KA_BETA_SPARSE_SIZE = 144`` + + KA_GAMMA_SPARSE_SIZE, : int + ``KA_GAMMA_SPARSE_SIZE = 144`` + + CA_SPARSE_SIZE, : int + ``CA_SPARSE_SIZE = 144`` + """ cdef public int KC0_SPARSE_SIZE cdef public int KG_SPARSE_SIZE cdef public int M_SPARSE_SIZE + cdef public int KA_BETA_SPARSE_SIZE + cdef public int KA_GAMMA_SPARSE_SIZE + cdef public int CA_SPARSE_SIZE def __cinit__(Quad4RData self): self.KC0_SPARSE_SIZE = 576 self.KG_SPARSE_SIZE = 144 self.M_SPARSE_SIZE = 480 + self.KA_BETA_SPARSE_SIZE = 144 + self.KA_GAMMA_SPARSE_SIZE = 144 + self.CA_SPARSE_SIZE = 144 + cdef class Quad4RProbe: r""" @@ -72,6 +88,7 @@ cdef class Quad4RProbe: self.xe = np.zeros(NUM_NODES*DOF//2, dtype=np.float64) self.ue = np.zeros(NUM_NODES*DOF, dtype=np.float64) + cdef class Quad4R: r""" Nodal connectivity for the plate element similar to Nastran's CQUAD4:: @@ -125,6 +142,9 @@ cdef class Quad4R: init_k_KC0, init_k_KG, init_k_M : int Position in the arrays storing the sparse data for the structural matrices. + init_k_KA_beta, init_k_KA_gamma, init_k_CA : int + Position in the arrays storing the sparse data for the aerodynamic + matrices based on the Piston theory. probe, : :class:`.Quad4RProbe` object Pointer to the probe. @@ -133,6 +153,7 @@ cdef class Quad4R: cdef public int n1, n2, n3, n4 cdef public int c1, c2, c3, c4 cdef public int init_k_KC0, init_k_KG, init_k_M + cdef public int init_k_KA_beta, init_k_KA_gamma, init_k_CA cdef public double area cdef public double alphat # drilling penalty factor for stiffness matrix, see Eq. 2.20 in F.M. Adam, A.E. Mohamed, A.E. Hassaballa, Degenerated Four Nodes Shell Element with Drilling Degree of Freedom, IOSR J. Eng. 3 (2013) 10–20. www.iosrjen.org (accessed April 20, 2020). cdef public double r11, r12, r13, r21, r22, r23, r31, r32, r33 @@ -155,6 +176,9 @@ cdef class Quad4R: # self.init_k_KCNL = 0 self.init_k_KG = 0 self.init_k_M = 0 + self.init_k_KA_beta = 0 + self.init_k_KA_gamma = 0 + self.init_k_CA = 0 self.area = 0 self.alphat = 1. # based on recommended value of reference F.M. Adam, A.E. Mohamed, A.E. Hassaballa self.r11 = self.r12 = self.r13 = 0. @@ -5294,14 +5318,9 @@ cdef class Quad4R: cdef double r11, r12, r13, r21, r22, r23, r31, r32, r33 cdef double x1, x2, x3, x4 cdef double y1, y2, y3, y4 - cdef double j11, j12, j21, j22 - cdef double N1x, N2x, N3x, N4x - cdef double N1y, N2y, N3y, N4y - cdef double cxx, cyy, cxy cdef double h11, h12, h13, h14, h22, h23, h24, h33, h34, h44, valH1 cdef double xi, eta, wij, detJ, N1, N2, N3, N4 cdef double points[2] - cdef double weights[2] with nogil: intrho = prop.intrho @@ -11659,3 +11678,2415 @@ cdef class Quad4R: k += 1 Mv[k] += h44*intrhoz2*r31**2 + h44*intrhoz2*r32**2 + + cpdef void update_KA_beta(Quad4R self, + long [::1] KA_betar, + long [::1] KA_betac, + double [::1] KA_betav, + ): + r"""Update sparse vectors for piston-theory aerodynamic matrix `KA_{\beta}` + + Properties + ---------- + KA_betar : np.array + Array to store row positions of sparse values + KA_betac : np.array + Array to store column positions of sparse values + KA_betav : np.array + Array to store sparse values + + """ + cdef int c1, c2, c3, c4, i, j, k + cdef double r11, r12, r13, r21, r22, r23, r31, r32, r33 + cdef double x1, x2, x3, x4 + cdef double y1, y2, y3, y4 + cdef double j11, j12, j21, j22 + cdef double N1, N2, N3, N4 + cdef double N1x, N2x, N3x, N4x + cdef double N1y, N2y, N3y, N4y + cdef double xi, eta, wij, detJ + cdef double points[2] + + with nogil: + + # NOTE ignoring z in local coordinates + x1 = self.probe.xe[0] + y1 = self.probe.xe[1] + # z1 = self.probe.xe[2] + x2 = self.probe.xe[3] + y2 = self.probe.xe[4] + # z2 = self.probe.xe[5] + x3 = self.probe.xe[6] + y3 = self.probe.xe[7] + # z3 = self.probe.xe[8] + x4 = self.probe.xe[9] + y4 = self.probe.xe[10] + # z4 = self.probe.xe[11] + + # Local to global transformation + r11 = self.r11 + r12 = self.r12 + r13 = self.r13 + r21 = self.r21 + r22 = self.r22 + r23 = self.r23 + r31 = self.r31 + r32 = self.r32 + r33 = self.r33 + + # positions the global matrices + c1 = self.c1 + c2 = self.c2 + c3 = self.c3 + c4 = self.c4 + + k = self.init_k_KA_beta + KA_betar[k] = 0+c1 + KA_betac[k] = 0+c1 + k += 1 + KA_betar[k] = 0+c1 + KA_betac[k] = 1+c1 + k += 1 + KA_betar[k] = 0+c1 + KA_betac[k] = 2+c1 + k += 1 + KA_betar[k] = 0+c1 + KA_betac[k] = 0+c2 + k += 1 + KA_betar[k] = 0+c1 + KA_betac[k] = 1+c2 + k += 1 + KA_betar[k] = 0+c1 + KA_betac[k] = 2+c2 + k += 1 + KA_betar[k] = 0+c1 + KA_betac[k] = 0+c3 + k += 1 + KA_betar[k] = 0+c1 + KA_betac[k] = 1+c3 + k += 1 + KA_betar[k] = 0+c1 + KA_betac[k] = 2+c3 + k += 1 + KA_betar[k] = 0+c1 + KA_betac[k] = 0+c4 + k += 1 + KA_betar[k] = 0+c1 + KA_betac[k] = 1+c4 + k += 1 + KA_betar[k] = 0+c1 + KA_betac[k] = 2+c4 + k += 1 + KA_betar[k] = 1+c1 + KA_betac[k] = 0+c1 + k += 1 + KA_betar[k] = 1+c1 + KA_betac[k] = 1+c1 + k += 1 + KA_betar[k] = 1+c1 + KA_betac[k] = 2+c1 + k += 1 + KA_betar[k] = 1+c1 + KA_betac[k] = 0+c2 + k += 1 + KA_betar[k] = 1+c1 + KA_betac[k] = 1+c2 + k += 1 + KA_betar[k] = 1+c1 + KA_betac[k] = 2+c2 + k += 1 + KA_betar[k] = 1+c1 + KA_betac[k] = 0+c3 + k += 1 + KA_betar[k] = 1+c1 + KA_betac[k] = 1+c3 + k += 1 + KA_betar[k] = 1+c1 + KA_betac[k] = 2+c3 + k += 1 + KA_betar[k] = 1+c1 + KA_betac[k] = 0+c4 + k += 1 + KA_betar[k] = 1+c1 + KA_betac[k] = 1+c4 + k += 1 + KA_betar[k] = 1+c1 + KA_betac[k] = 2+c4 + k += 1 + KA_betar[k] = 2+c1 + KA_betac[k] = 0+c1 + k += 1 + KA_betar[k] = 2+c1 + KA_betac[k] = 1+c1 + k += 1 + KA_betar[k] = 2+c1 + KA_betac[k] = 2+c1 + k += 1 + KA_betar[k] = 2+c1 + KA_betac[k] = 0+c2 + k += 1 + KA_betar[k] = 2+c1 + KA_betac[k] = 1+c2 + k += 1 + KA_betar[k] = 2+c1 + KA_betac[k] = 2+c2 + k += 1 + KA_betar[k] = 2+c1 + KA_betac[k] = 0+c3 + k += 1 + KA_betar[k] = 2+c1 + KA_betac[k] = 1+c3 + k += 1 + KA_betar[k] = 2+c1 + KA_betac[k] = 2+c3 + k += 1 + KA_betar[k] = 2+c1 + KA_betac[k] = 0+c4 + k += 1 + KA_betar[k] = 2+c1 + KA_betac[k] = 1+c4 + k += 1 + KA_betar[k] = 2+c1 + KA_betac[k] = 2+c4 + k += 1 + KA_betar[k] = 0+c2 + KA_betac[k] = 0+c1 + k += 1 + KA_betar[k] = 0+c2 + KA_betac[k] = 1+c1 + k += 1 + KA_betar[k] = 0+c2 + KA_betac[k] = 2+c1 + k += 1 + KA_betar[k] = 0+c2 + KA_betac[k] = 0+c2 + k += 1 + KA_betar[k] = 0+c2 + KA_betac[k] = 1+c2 + k += 1 + KA_betar[k] = 0+c2 + KA_betac[k] = 2+c2 + k += 1 + KA_betar[k] = 0+c2 + KA_betac[k] = 0+c3 + k += 1 + KA_betar[k] = 0+c2 + KA_betac[k] = 1+c3 + k += 1 + KA_betar[k] = 0+c2 + KA_betac[k] = 2+c3 + k += 1 + KA_betar[k] = 0+c2 + KA_betac[k] = 0+c4 + k += 1 + KA_betar[k] = 0+c2 + KA_betac[k] = 1+c4 + k += 1 + KA_betar[k] = 0+c2 + KA_betac[k] = 2+c4 + k += 1 + KA_betar[k] = 1+c2 + KA_betac[k] = 0+c1 + k += 1 + KA_betar[k] = 1+c2 + KA_betac[k] = 1+c1 + k += 1 + KA_betar[k] = 1+c2 + KA_betac[k] = 2+c1 + k += 1 + KA_betar[k] = 1+c2 + KA_betac[k] = 0+c2 + k += 1 + KA_betar[k] = 1+c2 + KA_betac[k] = 1+c2 + k += 1 + KA_betar[k] = 1+c2 + KA_betac[k] = 2+c2 + k += 1 + KA_betar[k] = 1+c2 + KA_betac[k] = 0+c3 + k += 1 + KA_betar[k] = 1+c2 + KA_betac[k] = 1+c3 + k += 1 + KA_betar[k] = 1+c2 + KA_betac[k] = 2+c3 + k += 1 + KA_betar[k] = 1+c2 + KA_betac[k] = 0+c4 + k += 1 + KA_betar[k] = 1+c2 + KA_betac[k] = 1+c4 + k += 1 + KA_betar[k] = 1+c2 + KA_betac[k] = 2+c4 + k += 1 + KA_betar[k] = 2+c2 + KA_betac[k] = 0+c1 + k += 1 + KA_betar[k] = 2+c2 + KA_betac[k] = 1+c1 + k += 1 + KA_betar[k] = 2+c2 + KA_betac[k] = 2+c1 + k += 1 + KA_betar[k] = 2+c2 + KA_betac[k] = 0+c2 + k += 1 + KA_betar[k] = 2+c2 + KA_betac[k] = 1+c2 + k += 1 + KA_betar[k] = 2+c2 + KA_betac[k] = 2+c2 + k += 1 + KA_betar[k] = 2+c2 + KA_betac[k] = 0+c3 + k += 1 + KA_betar[k] = 2+c2 + KA_betac[k] = 1+c3 + k += 1 + KA_betar[k] = 2+c2 + KA_betac[k] = 2+c3 + k += 1 + KA_betar[k] = 2+c2 + KA_betac[k] = 0+c4 + k += 1 + KA_betar[k] = 2+c2 + KA_betac[k] = 1+c4 + k += 1 + KA_betar[k] = 2+c2 + KA_betac[k] = 2+c4 + k += 1 + KA_betar[k] = 0+c3 + KA_betac[k] = 0+c1 + k += 1 + KA_betar[k] = 0+c3 + KA_betac[k] = 1+c1 + k += 1 + KA_betar[k] = 0+c3 + KA_betac[k] = 2+c1 + k += 1 + KA_betar[k] = 0+c3 + KA_betac[k] = 0+c2 + k += 1 + KA_betar[k] = 0+c3 + KA_betac[k] = 1+c2 + k += 1 + KA_betar[k] = 0+c3 + KA_betac[k] = 2+c2 + k += 1 + KA_betar[k] = 0+c3 + KA_betac[k] = 0+c3 + k += 1 + KA_betar[k] = 0+c3 + KA_betac[k] = 1+c3 + k += 1 + KA_betar[k] = 0+c3 + KA_betac[k] = 2+c3 + k += 1 + KA_betar[k] = 0+c3 + KA_betac[k] = 0+c4 + k += 1 + KA_betar[k] = 0+c3 + KA_betac[k] = 1+c4 + k += 1 + KA_betar[k] = 0+c3 + KA_betac[k] = 2+c4 + k += 1 + KA_betar[k] = 1+c3 + KA_betac[k] = 0+c1 + k += 1 + KA_betar[k] = 1+c3 + KA_betac[k] = 1+c1 + k += 1 + KA_betar[k] = 1+c3 + KA_betac[k] = 2+c1 + k += 1 + KA_betar[k] = 1+c3 + KA_betac[k] = 0+c2 + k += 1 + KA_betar[k] = 1+c3 + KA_betac[k] = 1+c2 + k += 1 + KA_betar[k] = 1+c3 + KA_betac[k] = 2+c2 + k += 1 + KA_betar[k] = 1+c3 + KA_betac[k] = 0+c3 + k += 1 + KA_betar[k] = 1+c3 + KA_betac[k] = 1+c3 + k += 1 + KA_betar[k] = 1+c3 + KA_betac[k] = 2+c3 + k += 1 + KA_betar[k] = 1+c3 + KA_betac[k] = 0+c4 + k += 1 + KA_betar[k] = 1+c3 + KA_betac[k] = 1+c4 + k += 1 + KA_betar[k] = 1+c3 + KA_betac[k] = 2+c4 + k += 1 + KA_betar[k] = 2+c3 + KA_betac[k] = 0+c1 + k += 1 + KA_betar[k] = 2+c3 + KA_betac[k] = 1+c1 + k += 1 + KA_betar[k] = 2+c3 + KA_betac[k] = 2+c1 + k += 1 + KA_betar[k] = 2+c3 + KA_betac[k] = 0+c2 + k += 1 + KA_betar[k] = 2+c3 + KA_betac[k] = 1+c2 + k += 1 + KA_betar[k] = 2+c3 + KA_betac[k] = 2+c2 + k += 1 + KA_betar[k] = 2+c3 + KA_betac[k] = 0+c3 + k += 1 + KA_betar[k] = 2+c3 + KA_betac[k] = 1+c3 + k += 1 + KA_betar[k] = 2+c3 + KA_betac[k] = 2+c3 + k += 1 + KA_betar[k] = 2+c3 + KA_betac[k] = 0+c4 + k += 1 + KA_betar[k] = 2+c3 + KA_betac[k] = 1+c4 + k += 1 + KA_betar[k] = 2+c3 + KA_betac[k] = 2+c4 + k += 1 + KA_betar[k] = 0+c4 + KA_betac[k] = 0+c1 + k += 1 + KA_betar[k] = 0+c4 + KA_betac[k] = 1+c1 + k += 1 + KA_betar[k] = 0+c4 + KA_betac[k] = 2+c1 + k += 1 + KA_betar[k] = 0+c4 + KA_betac[k] = 0+c2 + k += 1 + KA_betar[k] = 0+c4 + KA_betac[k] = 1+c2 + k += 1 + KA_betar[k] = 0+c4 + KA_betac[k] = 2+c2 + k += 1 + KA_betar[k] = 0+c4 + KA_betac[k] = 0+c3 + k += 1 + KA_betar[k] = 0+c4 + KA_betac[k] = 1+c3 + k += 1 + KA_betar[k] = 0+c4 + KA_betac[k] = 2+c3 + k += 1 + KA_betar[k] = 0+c4 + KA_betac[k] = 0+c4 + k += 1 + KA_betar[k] = 0+c4 + KA_betac[k] = 1+c4 + k += 1 + KA_betar[k] = 0+c4 + KA_betac[k] = 2+c4 + k += 1 + KA_betar[k] = 1+c4 + KA_betac[k] = 0+c1 + k += 1 + KA_betar[k] = 1+c4 + KA_betac[k] = 1+c1 + k += 1 + KA_betar[k] = 1+c4 + KA_betac[k] = 2+c1 + k += 1 + KA_betar[k] = 1+c4 + KA_betac[k] = 0+c2 + k += 1 + KA_betar[k] = 1+c4 + KA_betac[k] = 1+c2 + k += 1 + KA_betar[k] = 1+c4 + KA_betac[k] = 2+c2 + k += 1 + KA_betar[k] = 1+c4 + KA_betac[k] = 0+c3 + k += 1 + KA_betar[k] = 1+c4 + KA_betac[k] = 1+c3 + k += 1 + KA_betar[k] = 1+c4 + KA_betac[k] = 2+c3 + k += 1 + KA_betar[k] = 1+c4 + KA_betac[k] = 0+c4 + k += 1 + KA_betar[k] = 1+c4 + KA_betac[k] = 1+c4 + k += 1 + KA_betar[k] = 1+c4 + KA_betac[k] = 2+c4 + k += 1 + KA_betar[k] = 2+c4 + KA_betac[k] = 0+c1 + k += 1 + KA_betar[k] = 2+c4 + KA_betac[k] = 1+c1 + k += 1 + KA_betar[k] = 2+c4 + KA_betac[k] = 2+c1 + k += 1 + KA_betar[k] = 2+c4 + KA_betac[k] = 0+c2 + k += 1 + KA_betar[k] = 2+c4 + KA_betac[k] = 1+c2 + k += 1 + KA_betar[k] = 2+c4 + KA_betac[k] = 2+c2 + k += 1 + KA_betar[k] = 2+c4 + KA_betac[k] = 0+c3 + k += 1 + KA_betar[k] = 2+c4 + KA_betac[k] = 1+c3 + k += 1 + KA_betar[k] = 2+c4 + KA_betac[k] = 2+c3 + k += 1 + KA_betar[k] = 2+c4 + KA_betac[k] = 0+c4 + k += 1 + KA_betar[k] = 2+c4 + KA_betac[k] = 1+c4 + k += 1 + KA_betar[k] = 2+c4 + KA_betac[k] = 2+c4 + + # NOTE full integration for KG with two-point Gauss-Legendre quadrature + wij = 1. + points[0] = -0.5773502691896257645092 + points[1] = +0.5773502691896257645092 + + for i in range(2): + xi = points[i] + for j in range(2): + eta = points[j] + + detJ = -0.125*eta*x1*y2 + 0.125*eta*x1*y3 + 0.125*eta*x2*y1 - 0.125*eta*x2*y4 - 0.125*eta*x3*y1 + 0.125*eta*x3*y4 + 0.125*eta*x4*y2 - 0.125*eta*x4*y3 - 0.125*x1*xi*y3 + 0.125*x1*xi*y4 + 0.125*x1*y2 - 0.125*x1*y4 + 0.125*x2*xi*y3 - 0.125*x2*xi*y4 - 0.125*x2*y1 + 0.125*x2*y3 + 0.125*x3*xi*y1 - 0.125*x3*xi*y2 - 0.125*x3*y2 + 0.125*x3*y4 - 0.125*x4*xi*y1 + 0.125*x4*xi*y2 + 0.125*x4*y1 - 0.125*x4*y3 + + j11 = 2.0*(-xi*y1 + xi*y2 - xi*y3 + xi*y4 + y1 + y2 - y3 - y4)/(eta*x1*y2 - eta*x1*y3 - eta*x2*y1 + eta*x2*y4 + eta*x3*y1 - eta*x3*y4 - eta*x4*y2 + eta*x4*y3 + x1*xi*y3 - x1*xi*y4 - x1*y2 + x1*y4 - x2*xi*y3 + x2*xi*y4 + x2*y1 - x2*y3 - x3*xi*y1 + x3*xi*y2 + x3*y2 - x3*y4 + x4*xi*y1 - x4*xi*y2 - x4*y1 + x4*y3) + j12 = 2.0*(eta*y1 - eta*y2 + eta*y3 - eta*y4 - y1 + y2 + y3 - y4)/(eta*x1*y2 - eta*x1*y3 - eta*x2*y1 + eta*x2*y4 + eta*x3*y1 - eta*x3*y4 - eta*x4*y2 + eta*x4*y3 + x1*xi*y3 - x1*xi*y4 - x1*y2 + x1*y4 - x2*xi*y3 + x2*xi*y4 + x2*y1 - x2*y3 - x3*xi*y1 + x3*xi*y2 + x3*y2 - x3*y4 + x4*xi*y1 - x4*xi*y2 - x4*y1 + x4*y3) + j21 = 2.0*(x1*xi - x1 - x2*xi - x2 + x3*xi + x3 - x4*xi + x4)/(eta*x1*y2 - eta*x1*y3 - eta*x2*y1 + eta*x2*y4 + eta*x3*y1 - eta*x3*y4 - eta*x4*y2 + eta*x4*y3 + x1*xi*y3 - x1*xi*y4 - x1*y2 + x1*y4 - x2*xi*y3 + x2*xi*y4 + x2*y1 - x2*y3 - x3*xi*y1 + x3*xi*y2 + x3*y2 - x3*y4 + x4*xi*y1 - x4*xi*y2 - x4*y1 + x4*y3) + j22 = 2.0*(-eta*x1 + eta*x2 - eta*x3 + eta*x4 + x1 - x2 - x3 + x4)/(eta*x1*y2 - eta*x1*y3 - eta*x2*y1 + eta*x2*y4 + eta*x3*y1 - eta*x3*y4 - eta*x4*y2 + eta*x4*y3 + x1*xi*y3 - x1*xi*y4 - x1*y2 + x1*y4 - x2*xi*y3 + x2*xi*y4 + x2*y1 - x2*y3 - x3*xi*y1 + x3*xi*y2 + x3*y2 - x3*y4 + x4*xi*y1 - x4*xi*y2 - x4*y1 + x4*y3) + + N1 = 0.25*eta*xi - 0.25*eta - 0.25*xi + 0.25 + N2 = -0.25*eta*xi - 0.25*eta + 0.25*xi + 0.25 + N3 = 0.25*eta*xi + 0.25*eta + 0.25*xi + 0.25 + N4 = -0.25*eta*xi + 0.25*eta - 0.25*xi + 0.25 + + N1x = 0.25*j11*(eta - 1) + 0.25*j12*(xi - 1) + N2x = -0.25*eta*j11 + 0.25*j11 - 0.25*j12*xi - 0.25*j12 + N3x = 0.25*j11*(eta + 1) + 0.25*j12*(xi + 1) + N4x = -0.25*eta*j11 - 0.25*j11 - 0.25*j12*xi + 0.25*j12 + + N1y = 0.25*j21*(eta - 1) + 0.25*j22*(xi - 1) + N2y = -0.25*eta*j21 + 0.25*j21 - 0.25*j22*xi - 0.25*j22 + N3y = 0.25*j21*(eta + 1) + 0.25*j22*(xi + 1) + N4y = -0.25*eta*j21 - 0.25*j21 - 0.25*j22*xi + 0.25*j22 + + k = self.init_k_KA_beta + KA_betav[k] += -N1*detJ*r13**2*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r23*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13**2*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r23*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13**2*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r23*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13**2*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r23*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r23*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r23**2*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r23*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r23*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r23**2*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r23*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r23*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r23**2*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r23*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r23*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r23**2*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r23*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r23*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r33**2*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r23*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r33**2*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r23*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r33**2*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r13*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r23*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N1*detJ*r33**2*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13**2*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r23*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13**2*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r23*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13**2*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r23*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13**2*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r23*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r23*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r23**2*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r23*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r23*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r23**2*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r23*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r23*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r23**2*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r23*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r23*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r23**2*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r23*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r23*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r33**2*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r23*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r33**2*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r23*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r33**2*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r13*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r23*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N2*detJ*r33**2*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13**2*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r23*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13**2*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r23*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13**2*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r23*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13**2*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r23*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r23*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r23**2*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r23*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r23*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r23**2*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r23*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r23*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r23**2*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r23*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r23*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r23**2*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r23*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r23*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r33**2*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r23*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r33**2*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r23*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r33**2*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r13*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r23*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N3*detJ*r33**2*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13**2*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r23*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13**2*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r23*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13**2*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r23*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13**2*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r23*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r23*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r23**2*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r23*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r23*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r23**2*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r23*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r23*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r23**2*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r23*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r23*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r23**2*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r23*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r23*r33*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r33**2*wij*(N1x*r11 + N1y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r23*r33*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r33**2*wij*(N2x*r11 + N2y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r23*r33*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r33**2*wij*(N3x*r11 + N3y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r13*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r23*r33*wij*(N4x*r11 + N4y*r21) + k += 1 + KA_betav[k] += -N4*detJ*r33**2*wij*(N4x*r11 + N4y*r21) + + + cpdef void update_KA_gamma(Quad4R self, + long [::1] KA_gammar, + long [::1] KA_gammac, + double [::1] KA_gammav, + ): + r"""Update sparse vectors for piston-theory aerodynamic matrix `KA_{\gamma}` + + Properties + ---------- + KA_gammar : np.array + Array to store row positions of sparse values + KA_gammac : np.array + Array to store column positions of sparse values + KA_gammav : np.array + Array to store sparse values + + """ + cdef int c1, c2, c3, c4, i, j, k + cdef double r11, r12, r13, r21, r22, r23, r31, r32, r33 + cdef double x1, x2, x3, x4 + cdef double y1, y2, y3, y4 + cdef double N1, N2, N3, N4 + cdef double xi, eta, wij, detJ + cdef double points[2] + + with nogil: + + # NOTE ignoring z in local coordinates + x1 = self.probe.xe[0] + y1 = self.probe.xe[1] + # z1 = self.probe.xe[2] + x2 = self.probe.xe[3] + y2 = self.probe.xe[4] + # z2 = self.probe.xe[5] + x3 = self.probe.xe[6] + y3 = self.probe.xe[7] + # z3 = self.probe.xe[8] + x4 = self.probe.xe[9] + y4 = self.probe.xe[10] + # z4 = self.probe.xe[11] + + # Local to global transformation + r11 = self.r11 + r12 = self.r12 + r13 = self.r13 + r21 = self.r21 + r22 = self.r22 + r23 = self.r23 + r31 = self.r31 + r32 = self.r32 + r33 = self.r33 + + # positions the global matrices + c1 = self.c1 + c2 = self.c2 + c3 = self.c3 + c4 = self.c4 + + k = self.init_k_KA_gamma + KA_gammar[k] = 0+c1 + KA_gammac[k] = 0+c1 + k += 1 + KA_gammar[k] = 0+c1 + KA_gammac[k] = 1+c1 + k += 1 + KA_gammar[k] = 0+c1 + KA_gammac[k] = 2+c1 + k += 1 + KA_gammar[k] = 0+c1 + KA_gammac[k] = 0+c2 + k += 1 + KA_gammar[k] = 0+c1 + KA_gammac[k] = 1+c2 + k += 1 + KA_gammar[k] = 0+c1 + KA_gammac[k] = 2+c2 + k += 1 + KA_gammar[k] = 0+c1 + KA_gammac[k] = 0+c3 + k += 1 + KA_gammar[k] = 0+c1 + KA_gammac[k] = 1+c3 + k += 1 + KA_gammar[k] = 0+c1 + KA_gammac[k] = 2+c3 + k += 1 + KA_gammar[k] = 0+c1 + KA_gammac[k] = 0+c4 + k += 1 + KA_gammar[k] = 0+c1 + KA_gammac[k] = 1+c4 + k += 1 + KA_gammar[k] = 0+c1 + KA_gammac[k] = 2+c4 + k += 1 + KA_gammar[k] = 1+c1 + KA_gammac[k] = 0+c1 + k += 1 + KA_gammar[k] = 1+c1 + KA_gammac[k] = 1+c1 + k += 1 + KA_gammar[k] = 1+c1 + KA_gammac[k] = 2+c1 + k += 1 + KA_gammar[k] = 1+c1 + KA_gammac[k] = 0+c2 + k += 1 + KA_gammar[k] = 1+c1 + KA_gammac[k] = 1+c2 + k += 1 + KA_gammar[k] = 1+c1 + KA_gammac[k] = 2+c2 + k += 1 + KA_gammar[k] = 1+c1 + KA_gammac[k] = 0+c3 + k += 1 + KA_gammar[k] = 1+c1 + KA_gammac[k] = 1+c3 + k += 1 + KA_gammar[k] = 1+c1 + KA_gammac[k] = 2+c3 + k += 1 + KA_gammar[k] = 1+c1 + KA_gammac[k] = 0+c4 + k += 1 + KA_gammar[k] = 1+c1 + KA_gammac[k] = 1+c4 + k += 1 + KA_gammar[k] = 1+c1 + KA_gammac[k] = 2+c4 + k += 1 + KA_gammar[k] = 2+c1 + KA_gammac[k] = 0+c1 + k += 1 + KA_gammar[k] = 2+c1 + KA_gammac[k] = 1+c1 + k += 1 + KA_gammar[k] = 2+c1 + KA_gammac[k] = 2+c1 + k += 1 + KA_gammar[k] = 2+c1 + KA_gammac[k] = 0+c2 + k += 1 + KA_gammar[k] = 2+c1 + KA_gammac[k] = 1+c2 + k += 1 + KA_gammar[k] = 2+c1 + KA_gammac[k] = 2+c2 + k += 1 + KA_gammar[k] = 2+c1 + KA_gammac[k] = 0+c3 + k += 1 + KA_gammar[k] = 2+c1 + KA_gammac[k] = 1+c3 + k += 1 + KA_gammar[k] = 2+c1 + KA_gammac[k] = 2+c3 + k += 1 + KA_gammar[k] = 2+c1 + KA_gammac[k] = 0+c4 + k += 1 + KA_gammar[k] = 2+c1 + KA_gammac[k] = 1+c4 + k += 1 + KA_gammar[k] = 2+c1 + KA_gammac[k] = 2+c4 + k += 1 + KA_gammar[k] = 0+c2 + KA_gammac[k] = 0+c1 + k += 1 + KA_gammar[k] = 0+c2 + KA_gammac[k] = 1+c1 + k += 1 + KA_gammar[k] = 0+c2 + KA_gammac[k] = 2+c1 + k += 1 + KA_gammar[k] = 0+c2 + KA_gammac[k] = 0+c2 + k += 1 + KA_gammar[k] = 0+c2 + KA_gammac[k] = 1+c2 + k += 1 + KA_gammar[k] = 0+c2 + KA_gammac[k] = 2+c2 + k += 1 + KA_gammar[k] = 0+c2 + KA_gammac[k] = 0+c3 + k += 1 + KA_gammar[k] = 0+c2 + KA_gammac[k] = 1+c3 + k += 1 + KA_gammar[k] = 0+c2 + KA_gammac[k] = 2+c3 + k += 1 + KA_gammar[k] = 0+c2 + KA_gammac[k] = 0+c4 + k += 1 + KA_gammar[k] = 0+c2 + KA_gammac[k] = 1+c4 + k += 1 + KA_gammar[k] = 0+c2 + KA_gammac[k] = 2+c4 + k += 1 + KA_gammar[k] = 1+c2 + KA_gammac[k] = 0+c1 + k += 1 + KA_gammar[k] = 1+c2 + KA_gammac[k] = 1+c1 + k += 1 + KA_gammar[k] = 1+c2 + KA_gammac[k] = 2+c1 + k += 1 + KA_gammar[k] = 1+c2 + KA_gammac[k] = 0+c2 + k += 1 + KA_gammar[k] = 1+c2 + KA_gammac[k] = 1+c2 + k += 1 + KA_gammar[k] = 1+c2 + KA_gammac[k] = 2+c2 + k += 1 + KA_gammar[k] = 1+c2 + KA_gammac[k] = 0+c3 + k += 1 + KA_gammar[k] = 1+c2 + KA_gammac[k] = 1+c3 + k += 1 + KA_gammar[k] = 1+c2 + KA_gammac[k] = 2+c3 + k += 1 + KA_gammar[k] = 1+c2 + KA_gammac[k] = 0+c4 + k += 1 + KA_gammar[k] = 1+c2 + KA_gammac[k] = 1+c4 + k += 1 + KA_gammar[k] = 1+c2 + KA_gammac[k] = 2+c4 + k += 1 + KA_gammar[k] = 2+c2 + KA_gammac[k] = 0+c1 + k += 1 + KA_gammar[k] = 2+c2 + KA_gammac[k] = 1+c1 + k += 1 + KA_gammar[k] = 2+c2 + KA_gammac[k] = 2+c1 + k += 1 + KA_gammar[k] = 2+c2 + KA_gammac[k] = 0+c2 + k += 1 + KA_gammar[k] = 2+c2 + KA_gammac[k] = 1+c2 + k += 1 + KA_gammar[k] = 2+c2 + KA_gammac[k] = 2+c2 + k += 1 + KA_gammar[k] = 2+c2 + KA_gammac[k] = 0+c3 + k += 1 + KA_gammar[k] = 2+c2 + KA_gammac[k] = 1+c3 + k += 1 + KA_gammar[k] = 2+c2 + KA_gammac[k] = 2+c3 + k += 1 + KA_gammar[k] = 2+c2 + KA_gammac[k] = 0+c4 + k += 1 + KA_gammar[k] = 2+c2 + KA_gammac[k] = 1+c4 + k += 1 + KA_gammar[k] = 2+c2 + KA_gammac[k] = 2+c4 + k += 1 + KA_gammar[k] = 0+c3 + KA_gammac[k] = 0+c1 + k += 1 + KA_gammar[k] = 0+c3 + KA_gammac[k] = 1+c1 + k += 1 + KA_gammar[k] = 0+c3 + KA_gammac[k] = 2+c1 + k += 1 + KA_gammar[k] = 0+c3 + KA_gammac[k] = 0+c2 + k += 1 + KA_gammar[k] = 0+c3 + KA_gammac[k] = 1+c2 + k += 1 + KA_gammar[k] = 0+c3 + KA_gammac[k] = 2+c2 + k += 1 + KA_gammar[k] = 0+c3 + KA_gammac[k] = 0+c3 + k += 1 + KA_gammar[k] = 0+c3 + KA_gammac[k] = 1+c3 + k += 1 + KA_gammar[k] = 0+c3 + KA_gammac[k] = 2+c3 + k += 1 + KA_gammar[k] = 0+c3 + KA_gammac[k] = 0+c4 + k += 1 + KA_gammar[k] = 0+c3 + KA_gammac[k] = 1+c4 + k += 1 + KA_gammar[k] = 0+c3 + KA_gammac[k] = 2+c4 + k += 1 + KA_gammar[k] = 1+c3 + KA_gammac[k] = 0+c1 + k += 1 + KA_gammar[k] = 1+c3 + KA_gammac[k] = 1+c1 + k += 1 + KA_gammar[k] = 1+c3 + KA_gammac[k] = 2+c1 + k += 1 + KA_gammar[k] = 1+c3 + KA_gammac[k] = 0+c2 + k += 1 + KA_gammar[k] = 1+c3 + KA_gammac[k] = 1+c2 + k += 1 + KA_gammar[k] = 1+c3 + KA_gammac[k] = 2+c2 + k += 1 + KA_gammar[k] = 1+c3 + KA_gammac[k] = 0+c3 + k += 1 + KA_gammar[k] = 1+c3 + KA_gammac[k] = 1+c3 + k += 1 + KA_gammar[k] = 1+c3 + KA_gammac[k] = 2+c3 + k += 1 + KA_gammar[k] = 1+c3 + KA_gammac[k] = 0+c4 + k += 1 + KA_gammar[k] = 1+c3 + KA_gammac[k] = 1+c4 + k += 1 + KA_gammar[k] = 1+c3 + KA_gammac[k] = 2+c4 + k += 1 + KA_gammar[k] = 2+c3 + KA_gammac[k] = 0+c1 + k += 1 + KA_gammar[k] = 2+c3 + KA_gammac[k] = 1+c1 + k += 1 + KA_gammar[k] = 2+c3 + KA_gammac[k] = 2+c1 + k += 1 + KA_gammar[k] = 2+c3 + KA_gammac[k] = 0+c2 + k += 1 + KA_gammar[k] = 2+c3 + KA_gammac[k] = 1+c2 + k += 1 + KA_gammar[k] = 2+c3 + KA_gammac[k] = 2+c2 + k += 1 + KA_gammar[k] = 2+c3 + KA_gammac[k] = 0+c3 + k += 1 + KA_gammar[k] = 2+c3 + KA_gammac[k] = 1+c3 + k += 1 + KA_gammar[k] = 2+c3 + KA_gammac[k] = 2+c3 + k += 1 + KA_gammar[k] = 2+c3 + KA_gammac[k] = 0+c4 + k += 1 + KA_gammar[k] = 2+c3 + KA_gammac[k] = 1+c4 + k += 1 + KA_gammar[k] = 2+c3 + KA_gammac[k] = 2+c4 + k += 1 + KA_gammar[k] = 0+c4 + KA_gammac[k] = 0+c1 + k += 1 + KA_gammar[k] = 0+c4 + KA_gammac[k] = 1+c1 + k += 1 + KA_gammar[k] = 0+c4 + KA_gammac[k] = 2+c1 + k += 1 + KA_gammar[k] = 0+c4 + KA_gammac[k] = 0+c2 + k += 1 + KA_gammar[k] = 0+c4 + KA_gammac[k] = 1+c2 + k += 1 + KA_gammar[k] = 0+c4 + KA_gammac[k] = 2+c2 + k += 1 + KA_gammar[k] = 0+c4 + KA_gammac[k] = 0+c3 + k += 1 + KA_gammar[k] = 0+c4 + KA_gammac[k] = 1+c3 + k += 1 + KA_gammar[k] = 0+c4 + KA_gammac[k] = 2+c3 + k += 1 + KA_gammar[k] = 0+c4 + KA_gammac[k] = 0+c4 + k += 1 + KA_gammar[k] = 0+c4 + KA_gammac[k] = 1+c4 + k += 1 + KA_gammar[k] = 0+c4 + KA_gammac[k] = 2+c4 + k += 1 + KA_gammar[k] = 1+c4 + KA_gammac[k] = 0+c1 + k += 1 + KA_gammar[k] = 1+c4 + KA_gammac[k] = 1+c1 + k += 1 + KA_gammar[k] = 1+c4 + KA_gammac[k] = 2+c1 + k += 1 + KA_gammar[k] = 1+c4 + KA_gammac[k] = 0+c2 + k += 1 + KA_gammar[k] = 1+c4 + KA_gammac[k] = 1+c2 + k += 1 + KA_gammar[k] = 1+c4 + KA_gammac[k] = 2+c2 + k += 1 + KA_gammar[k] = 1+c4 + KA_gammac[k] = 0+c3 + k += 1 + KA_gammar[k] = 1+c4 + KA_gammac[k] = 1+c3 + k += 1 + KA_gammar[k] = 1+c4 + KA_gammac[k] = 2+c3 + k += 1 + KA_gammar[k] = 1+c4 + KA_gammac[k] = 0+c4 + k += 1 + KA_gammar[k] = 1+c4 + KA_gammac[k] = 1+c4 + k += 1 + KA_gammar[k] = 1+c4 + KA_gammac[k] = 2+c4 + k += 1 + KA_gammar[k] = 2+c4 + KA_gammac[k] = 0+c1 + k += 1 + KA_gammar[k] = 2+c4 + KA_gammac[k] = 1+c1 + k += 1 + KA_gammar[k] = 2+c4 + KA_gammac[k] = 2+c1 + k += 1 + KA_gammar[k] = 2+c4 + KA_gammac[k] = 0+c2 + k += 1 + KA_gammar[k] = 2+c4 + KA_gammac[k] = 1+c2 + k += 1 + KA_gammar[k] = 2+c4 + KA_gammac[k] = 2+c2 + k += 1 + KA_gammar[k] = 2+c4 + KA_gammac[k] = 0+c3 + k += 1 + KA_gammar[k] = 2+c4 + KA_gammac[k] = 1+c3 + k += 1 + KA_gammar[k] = 2+c4 + KA_gammac[k] = 2+c3 + k += 1 + KA_gammar[k] = 2+c4 + KA_gammac[k] = 0+c4 + k += 1 + KA_gammar[k] = 2+c4 + KA_gammac[k] = 1+c4 + k += 1 + KA_gammar[k] = 2+c4 + KA_gammac[k] = 2+c4 + + # NOTE full integration for KG with two-point Gauss-Legendre quadrature + wij = 1. + points[0] = -0.5773502691896257645092 + points[1] = +0.5773502691896257645092 + + for i in range(2): + xi = points[i] + for j in range(2): + eta = points[j] + + detJ = -0.125*eta*x1*y2 + 0.125*eta*x1*y3 + 0.125*eta*x2*y1 - 0.125*eta*x2*y4 - 0.125*eta*x3*y1 + 0.125*eta*x3*y4 + 0.125*eta*x4*y2 - 0.125*eta*x4*y3 - 0.125*x1*xi*y3 + 0.125*x1*xi*y4 + 0.125*x1*y2 - 0.125*x1*y4 + 0.125*x2*xi*y3 - 0.125*x2*xi*y4 - 0.125*x2*y1 + 0.125*x2*y3 + 0.125*x3*xi*y1 - 0.125*x3*xi*y2 - 0.125*x3*y2 + 0.125*x3*y4 - 0.125*x4*xi*y1 + 0.125*x4*xi*y2 + 0.125*x4*y1 - 0.125*x4*y3 + + N1 = 0.25*eta*xi - 0.25*eta - 0.25*xi + 0.25 + N2 = -0.25*eta*xi - 0.25*eta + 0.25*xi + 0.25 + N3 = 0.25*eta*xi + 0.25*eta + 0.25*xi + 0.25 + N4 = -0.25*eta*xi + 0.25*eta - 0.25*xi + 0.25 + + k = self.init_k_KA_gamma + KA_gammav[k] += N1**2*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N1**2*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N1**2*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N1**2*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N1**2*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N1**2*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N1**2*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N1**2*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N1**2*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N2**2*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N2**2*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N2**2*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N2**2*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N2**2*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N2**2*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N1*N2*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N2**2*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N2**2*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N2**2*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N3**2*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N3**2*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N3**2*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N3**2*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N3**2*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N3**2*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N1*N3*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N2*N3*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N3**2*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N3**2*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N3**2*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N4**2*detJ*r13**2*wij + k += 1 + KA_gammav[k] += N4**2*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N4**2*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N4**2*detJ*r13*r23*wij + k += 1 + KA_gammav[k] += N4**2*detJ*r23**2*wij + k += 1 + KA_gammav[k] += N4**2*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N1*N4*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N2*N4*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N3*N4*detJ*r33**2*wij + k += 1 + KA_gammav[k] += N4**2*detJ*r13*r33*wij + k += 1 + KA_gammav[k] += N4**2*detJ*r23*r33*wij + k += 1 + KA_gammav[k] += N4**2*detJ*r33**2*wij + + + cpdef void update_CA(Quad4R self, + long [::1] CAr, + long [::1] CAc, + double [::1] CAv, + ): + r"""Update sparse vectors for piston-theory aerodynamic damping matrix `CA` + + Properties + ---------- + CAr : np.array + Array to store row positions of sparse values + CAc : np.array + Array to store column positions of sparse values + CAv : np.array + Array to store sparse values + + """ + cdef int c1, c2, c3, c4, i, j, k + cdef double r11, r12, r13, r21, r22, r23, r31, r32, r33 + cdef double x1, x2, x3, x4 + cdef double y1, y2, y3, y4 + cdef double N1, N2, N3, N4 + cdef double xi, eta, wij, detJ + cdef double points[2] + + with nogil: + + # NOTE ignoring z in local coordinates + x1 = self.probe.xe[0] + y1 = self.probe.xe[1] + # z1 = self.probe.xe[2] + x2 = self.probe.xe[3] + y2 = self.probe.xe[4] + # z2 = self.probe.xe[5] + x3 = self.probe.xe[6] + y3 = self.probe.xe[7] + # z3 = self.probe.xe[8] + x4 = self.probe.xe[9] + y4 = self.probe.xe[10] + # z4 = self.probe.xe[11] + + # Local to global transformation + r11 = self.r11 + r12 = self.r12 + r13 = self.r13 + r21 = self.r21 + r22 = self.r22 + r23 = self.r23 + r31 = self.r31 + r32 = self.r32 + r33 = self.r33 + + # positions the global matrices + c1 = self.c1 + c2 = self.c2 + c3 = self.c3 + c4 = self.c4 + + k = self.init_k_CA + CAr[k] = 0+c1 + CAc[k] = 0+c1 + k += 1 + CAr[k] = 0+c1 + CAc[k] = 1+c1 + k += 1 + CAr[k] = 0+c1 + CAc[k] = 2+c1 + k += 1 + CAr[k] = 0+c1 + CAc[k] = 0+c2 + k += 1 + CAr[k] = 0+c1 + CAc[k] = 1+c2 + k += 1 + CAr[k] = 0+c1 + CAc[k] = 2+c2 + k += 1 + CAr[k] = 0+c1 + CAc[k] = 0+c3 + k += 1 + CAr[k] = 0+c1 + CAc[k] = 1+c3 + k += 1 + CAr[k] = 0+c1 + CAc[k] = 2+c3 + k += 1 + CAr[k] = 0+c1 + CAc[k] = 0+c4 + k += 1 + CAr[k] = 0+c1 + CAc[k] = 1+c4 + k += 1 + CAr[k] = 0+c1 + CAc[k] = 2+c4 + k += 1 + CAr[k] = 1+c1 + CAc[k] = 0+c1 + k += 1 + CAr[k] = 1+c1 + CAc[k] = 1+c1 + k += 1 + CAr[k] = 1+c1 + CAc[k] = 2+c1 + k += 1 + CAr[k] = 1+c1 + CAc[k] = 0+c2 + k += 1 + CAr[k] = 1+c1 + CAc[k] = 1+c2 + k += 1 + CAr[k] = 1+c1 + CAc[k] = 2+c2 + k += 1 + CAr[k] = 1+c1 + CAc[k] = 0+c3 + k += 1 + CAr[k] = 1+c1 + CAc[k] = 1+c3 + k += 1 + CAr[k] = 1+c1 + CAc[k] = 2+c3 + k += 1 + CAr[k] = 1+c1 + CAc[k] = 0+c4 + k += 1 + CAr[k] = 1+c1 + CAc[k] = 1+c4 + k += 1 + CAr[k] = 1+c1 + CAc[k] = 2+c4 + k += 1 + CAr[k] = 2+c1 + CAc[k] = 0+c1 + k += 1 + CAr[k] = 2+c1 + CAc[k] = 1+c1 + k += 1 + CAr[k] = 2+c1 + CAc[k] = 2+c1 + k += 1 + CAr[k] = 2+c1 + CAc[k] = 0+c2 + k += 1 + CAr[k] = 2+c1 + CAc[k] = 1+c2 + k += 1 + CAr[k] = 2+c1 + CAc[k] = 2+c2 + k += 1 + CAr[k] = 2+c1 + CAc[k] = 0+c3 + k += 1 + CAr[k] = 2+c1 + CAc[k] = 1+c3 + k += 1 + CAr[k] = 2+c1 + CAc[k] = 2+c3 + k += 1 + CAr[k] = 2+c1 + CAc[k] = 0+c4 + k += 1 + CAr[k] = 2+c1 + CAc[k] = 1+c4 + k += 1 + CAr[k] = 2+c1 + CAc[k] = 2+c4 + k += 1 + CAr[k] = 0+c2 + CAc[k] = 0+c1 + k += 1 + CAr[k] = 0+c2 + CAc[k] = 1+c1 + k += 1 + CAr[k] = 0+c2 + CAc[k] = 2+c1 + k += 1 + CAr[k] = 0+c2 + CAc[k] = 0+c2 + k += 1 + CAr[k] = 0+c2 + CAc[k] = 1+c2 + k += 1 + CAr[k] = 0+c2 + CAc[k] = 2+c2 + k += 1 + CAr[k] = 0+c2 + CAc[k] = 0+c3 + k += 1 + CAr[k] = 0+c2 + CAc[k] = 1+c3 + k += 1 + CAr[k] = 0+c2 + CAc[k] = 2+c3 + k += 1 + CAr[k] = 0+c2 + CAc[k] = 0+c4 + k += 1 + CAr[k] = 0+c2 + CAc[k] = 1+c4 + k += 1 + CAr[k] = 0+c2 + CAc[k] = 2+c4 + k += 1 + CAr[k] = 1+c2 + CAc[k] = 0+c1 + k += 1 + CAr[k] = 1+c2 + CAc[k] = 1+c1 + k += 1 + CAr[k] = 1+c2 + CAc[k] = 2+c1 + k += 1 + CAr[k] = 1+c2 + CAc[k] = 0+c2 + k += 1 + CAr[k] = 1+c2 + CAc[k] = 1+c2 + k += 1 + CAr[k] = 1+c2 + CAc[k] = 2+c2 + k += 1 + CAr[k] = 1+c2 + CAc[k] = 0+c3 + k += 1 + CAr[k] = 1+c2 + CAc[k] = 1+c3 + k += 1 + CAr[k] = 1+c2 + CAc[k] = 2+c3 + k += 1 + CAr[k] = 1+c2 + CAc[k] = 0+c4 + k += 1 + CAr[k] = 1+c2 + CAc[k] = 1+c4 + k += 1 + CAr[k] = 1+c2 + CAc[k] = 2+c4 + k += 1 + CAr[k] = 2+c2 + CAc[k] = 0+c1 + k += 1 + CAr[k] = 2+c2 + CAc[k] = 1+c1 + k += 1 + CAr[k] = 2+c2 + CAc[k] = 2+c1 + k += 1 + CAr[k] = 2+c2 + CAc[k] = 0+c2 + k += 1 + CAr[k] = 2+c2 + CAc[k] = 1+c2 + k += 1 + CAr[k] = 2+c2 + CAc[k] = 2+c2 + k += 1 + CAr[k] = 2+c2 + CAc[k] = 0+c3 + k += 1 + CAr[k] = 2+c2 + CAc[k] = 1+c3 + k += 1 + CAr[k] = 2+c2 + CAc[k] = 2+c3 + k += 1 + CAr[k] = 2+c2 + CAc[k] = 0+c4 + k += 1 + CAr[k] = 2+c2 + CAc[k] = 1+c4 + k += 1 + CAr[k] = 2+c2 + CAc[k] = 2+c4 + k += 1 + CAr[k] = 0+c3 + CAc[k] = 0+c1 + k += 1 + CAr[k] = 0+c3 + CAc[k] = 1+c1 + k += 1 + CAr[k] = 0+c3 + CAc[k] = 2+c1 + k += 1 + CAr[k] = 0+c3 + CAc[k] = 0+c2 + k += 1 + CAr[k] = 0+c3 + CAc[k] = 1+c2 + k += 1 + CAr[k] = 0+c3 + CAc[k] = 2+c2 + k += 1 + CAr[k] = 0+c3 + CAc[k] = 0+c3 + k += 1 + CAr[k] = 0+c3 + CAc[k] = 1+c3 + k += 1 + CAr[k] = 0+c3 + CAc[k] = 2+c3 + k += 1 + CAr[k] = 0+c3 + CAc[k] = 0+c4 + k += 1 + CAr[k] = 0+c3 + CAc[k] = 1+c4 + k += 1 + CAr[k] = 0+c3 + CAc[k] = 2+c4 + k += 1 + CAr[k] = 1+c3 + CAc[k] = 0+c1 + k += 1 + CAr[k] = 1+c3 + CAc[k] = 1+c1 + k += 1 + CAr[k] = 1+c3 + CAc[k] = 2+c1 + k += 1 + CAr[k] = 1+c3 + CAc[k] = 0+c2 + k += 1 + CAr[k] = 1+c3 + CAc[k] = 1+c2 + k += 1 + CAr[k] = 1+c3 + CAc[k] = 2+c2 + k += 1 + CAr[k] = 1+c3 + CAc[k] = 0+c3 + k += 1 + CAr[k] = 1+c3 + CAc[k] = 1+c3 + k += 1 + CAr[k] = 1+c3 + CAc[k] = 2+c3 + k += 1 + CAr[k] = 1+c3 + CAc[k] = 0+c4 + k += 1 + CAr[k] = 1+c3 + CAc[k] = 1+c4 + k += 1 + CAr[k] = 1+c3 + CAc[k] = 2+c4 + k += 1 + CAr[k] = 2+c3 + CAc[k] = 0+c1 + k += 1 + CAr[k] = 2+c3 + CAc[k] = 1+c1 + k += 1 + CAr[k] = 2+c3 + CAc[k] = 2+c1 + k += 1 + CAr[k] = 2+c3 + CAc[k] = 0+c2 + k += 1 + CAr[k] = 2+c3 + CAc[k] = 1+c2 + k += 1 + CAr[k] = 2+c3 + CAc[k] = 2+c2 + k += 1 + CAr[k] = 2+c3 + CAc[k] = 0+c3 + k += 1 + CAr[k] = 2+c3 + CAc[k] = 1+c3 + k += 1 + CAr[k] = 2+c3 + CAc[k] = 2+c3 + k += 1 + CAr[k] = 2+c3 + CAc[k] = 0+c4 + k += 1 + CAr[k] = 2+c3 + CAc[k] = 1+c4 + k += 1 + CAr[k] = 2+c3 + CAc[k] = 2+c4 + k += 1 + CAr[k] = 0+c4 + CAc[k] = 0+c1 + k += 1 + CAr[k] = 0+c4 + CAc[k] = 1+c1 + k += 1 + CAr[k] = 0+c4 + CAc[k] = 2+c1 + k += 1 + CAr[k] = 0+c4 + CAc[k] = 0+c2 + k += 1 + CAr[k] = 0+c4 + CAc[k] = 1+c2 + k += 1 + CAr[k] = 0+c4 + CAc[k] = 2+c2 + k += 1 + CAr[k] = 0+c4 + CAc[k] = 0+c3 + k += 1 + CAr[k] = 0+c4 + CAc[k] = 1+c3 + k += 1 + CAr[k] = 0+c4 + CAc[k] = 2+c3 + k += 1 + CAr[k] = 0+c4 + CAc[k] = 0+c4 + k += 1 + CAr[k] = 0+c4 + CAc[k] = 1+c4 + k += 1 + CAr[k] = 0+c4 + CAc[k] = 2+c4 + k += 1 + CAr[k] = 1+c4 + CAc[k] = 0+c1 + k += 1 + CAr[k] = 1+c4 + CAc[k] = 1+c1 + k += 1 + CAr[k] = 1+c4 + CAc[k] = 2+c1 + k += 1 + CAr[k] = 1+c4 + CAc[k] = 0+c2 + k += 1 + CAr[k] = 1+c4 + CAc[k] = 1+c2 + k += 1 + CAr[k] = 1+c4 + CAc[k] = 2+c2 + k += 1 + CAr[k] = 1+c4 + CAc[k] = 0+c3 + k += 1 + CAr[k] = 1+c4 + CAc[k] = 1+c3 + k += 1 + CAr[k] = 1+c4 + CAc[k] = 2+c3 + k += 1 + CAr[k] = 1+c4 + CAc[k] = 0+c4 + k += 1 + CAr[k] = 1+c4 + CAc[k] = 1+c4 + k += 1 + CAr[k] = 1+c4 + CAc[k] = 2+c4 + k += 1 + CAr[k] = 2+c4 + CAc[k] = 0+c1 + k += 1 + CAr[k] = 2+c4 + CAc[k] = 1+c1 + k += 1 + CAr[k] = 2+c4 + CAc[k] = 2+c1 + k += 1 + CAr[k] = 2+c4 + CAc[k] = 0+c2 + k += 1 + CAr[k] = 2+c4 + CAc[k] = 1+c2 + k += 1 + CAr[k] = 2+c4 + CAc[k] = 2+c2 + k += 1 + CAr[k] = 2+c4 + CAc[k] = 0+c3 + k += 1 + CAr[k] = 2+c4 + CAc[k] = 1+c3 + k += 1 + CAr[k] = 2+c4 + CAc[k] = 2+c3 + k += 1 + CAr[k] = 2+c4 + CAc[k] = 0+c4 + k += 1 + CAr[k] = 2+c4 + CAc[k] = 1+c4 + k += 1 + CAr[k] = 2+c4 + CAc[k] = 2+c4 + + # NOTE full integration for KG with two-point Gauss-Legendre quadrature + wij = 1. + points[0] = -0.5773502691896257645092 + points[1] = +0.5773502691896257645092 + + for i in range(2): + xi = points[i] + for j in range(2): + eta = points[j] + + detJ = -0.125*eta*x1*y2 + 0.125*eta*x1*y3 + 0.125*eta*x2*y1 - 0.125*eta*x2*y4 - 0.125*eta*x3*y1 + 0.125*eta*x3*y4 + 0.125*eta*x4*y2 - 0.125*eta*x4*y3 - 0.125*x1*xi*y3 + 0.125*x1*xi*y4 + 0.125*x1*y2 - 0.125*x1*y4 + 0.125*x2*xi*y3 - 0.125*x2*xi*y4 - 0.125*x2*y1 + 0.125*x2*y3 + 0.125*x3*xi*y1 - 0.125*x3*xi*y2 - 0.125*x3*y2 + 0.125*x3*y4 - 0.125*x4*xi*y1 + 0.125*x4*xi*y2 + 0.125*x4*y1 - 0.125*x4*y3 + + N1 = 0.25*eta*xi - 0.25*eta - 0.25*xi + 0.25 + N2 = -0.25*eta*xi - 0.25*eta + 0.25*xi + 0.25 + N3 = 0.25*eta*xi + 0.25*eta + 0.25*xi + 0.25 + N4 = -0.25*eta*xi + 0.25*eta - 0.25*xi + 0.25 + + k = self.init_k_CA + CAv[k] += -N1**2*detJ*r13**2*wij + k += 1 + CAv[k] += -N1**2*detJ*r13*r23*wij + k += 1 + CAv[k] += -N1**2*detJ*r13*r33*wij + k += 1 + CAv[k] += -N1*N2*detJ*r13**2*wij + k += 1 + CAv[k] += -N1*N2*detJ*r13*r23*wij + k += 1 + CAv[k] += -N1*N2*detJ*r13*r33*wij + k += 1 + CAv[k] += -N1*N3*detJ*r13**2*wij + k += 1 + CAv[k] += -N1*N3*detJ*r13*r23*wij + k += 1 + CAv[k] += -N1*N3*detJ*r13*r33*wij + k += 1 + CAv[k] += -N1*N4*detJ*r13**2*wij + k += 1 + CAv[k] += -N1*N4*detJ*r13*r23*wij + k += 1 + CAv[k] += -N1*N4*detJ*r13*r33*wij + k += 1 + CAv[k] += -N1**2*detJ*r13*r23*wij + k += 1 + CAv[k] += -N1**2*detJ*r23**2*wij + k += 1 + CAv[k] += -N1**2*detJ*r23*r33*wij + k += 1 + CAv[k] += -N1*N2*detJ*r13*r23*wij + k += 1 + CAv[k] += -N1*N2*detJ*r23**2*wij + k += 1 + CAv[k] += -N1*N2*detJ*r23*r33*wij + k += 1 + CAv[k] += -N1*N3*detJ*r13*r23*wij + k += 1 + CAv[k] += -N1*N3*detJ*r23**2*wij + k += 1 + CAv[k] += -N1*N3*detJ*r23*r33*wij + k += 1 + CAv[k] += -N1*N4*detJ*r13*r23*wij + k += 1 + CAv[k] += -N1*N4*detJ*r23**2*wij + k += 1 + CAv[k] += -N1*N4*detJ*r23*r33*wij + k += 1 + CAv[k] += -N1**2*detJ*r13*r33*wij + k += 1 + CAv[k] += -N1**2*detJ*r23*r33*wij + k += 1 + CAv[k] += -N1**2*detJ*r33**2*wij + k += 1 + CAv[k] += -N1*N2*detJ*r13*r33*wij + k += 1 + CAv[k] += -N1*N2*detJ*r23*r33*wij + k += 1 + CAv[k] += -N1*N2*detJ*r33**2*wij + k += 1 + CAv[k] += -N1*N3*detJ*r13*r33*wij + k += 1 + CAv[k] += -N1*N3*detJ*r23*r33*wij + k += 1 + CAv[k] += -N1*N3*detJ*r33**2*wij + k += 1 + CAv[k] += -N1*N4*detJ*r13*r33*wij + k += 1 + CAv[k] += -N1*N4*detJ*r23*r33*wij + k += 1 + CAv[k] += -N1*N4*detJ*r33**2*wij + k += 1 + CAv[k] += -N1*N2*detJ*r13**2*wij + k += 1 + CAv[k] += -N1*N2*detJ*r13*r23*wij + k += 1 + CAv[k] += -N1*N2*detJ*r13*r33*wij + k += 1 + CAv[k] += -N2**2*detJ*r13**2*wij + k += 1 + CAv[k] += -N2**2*detJ*r13*r23*wij + k += 1 + CAv[k] += -N2**2*detJ*r13*r33*wij + k += 1 + CAv[k] += -N2*N3*detJ*r13**2*wij + k += 1 + CAv[k] += -N2*N3*detJ*r13*r23*wij + k += 1 + CAv[k] += -N2*N3*detJ*r13*r33*wij + k += 1 + CAv[k] += -N2*N4*detJ*r13**2*wij + k += 1 + CAv[k] += -N2*N4*detJ*r13*r23*wij + k += 1 + CAv[k] += -N2*N4*detJ*r13*r33*wij + k += 1 + CAv[k] += -N1*N2*detJ*r13*r23*wij + k += 1 + CAv[k] += -N1*N2*detJ*r23**2*wij + k += 1 + CAv[k] += -N1*N2*detJ*r23*r33*wij + k += 1 + CAv[k] += -N2**2*detJ*r13*r23*wij + k += 1 + CAv[k] += -N2**2*detJ*r23**2*wij + k += 1 + CAv[k] += -N2**2*detJ*r23*r33*wij + k += 1 + CAv[k] += -N2*N3*detJ*r13*r23*wij + k += 1 + CAv[k] += -N2*N3*detJ*r23**2*wij + k += 1 + CAv[k] += -N2*N3*detJ*r23*r33*wij + k += 1 + CAv[k] += -N2*N4*detJ*r13*r23*wij + k += 1 + CAv[k] += -N2*N4*detJ*r23**2*wij + k += 1 + CAv[k] += -N2*N4*detJ*r23*r33*wij + k += 1 + CAv[k] += -N1*N2*detJ*r13*r33*wij + k += 1 + CAv[k] += -N1*N2*detJ*r23*r33*wij + k += 1 + CAv[k] += -N1*N2*detJ*r33**2*wij + k += 1 + CAv[k] += -N2**2*detJ*r13*r33*wij + k += 1 + CAv[k] += -N2**2*detJ*r23*r33*wij + k += 1 + CAv[k] += -N2**2*detJ*r33**2*wij + k += 1 + CAv[k] += -N2*N3*detJ*r13*r33*wij + k += 1 + CAv[k] += -N2*N3*detJ*r23*r33*wij + k += 1 + CAv[k] += -N2*N3*detJ*r33**2*wij + k += 1 + CAv[k] += -N2*N4*detJ*r13*r33*wij + k += 1 + CAv[k] += -N2*N4*detJ*r23*r33*wij + k += 1 + CAv[k] += -N2*N4*detJ*r33**2*wij + k += 1 + CAv[k] += -N1*N3*detJ*r13**2*wij + k += 1 + CAv[k] += -N1*N3*detJ*r13*r23*wij + k += 1 + CAv[k] += -N1*N3*detJ*r13*r33*wij + k += 1 + CAv[k] += -N2*N3*detJ*r13**2*wij + k += 1 + CAv[k] += -N2*N3*detJ*r13*r23*wij + k += 1 + CAv[k] += -N2*N3*detJ*r13*r33*wij + k += 1 + CAv[k] += -N3**2*detJ*r13**2*wij + k += 1 + CAv[k] += -N3**2*detJ*r13*r23*wij + k += 1 + CAv[k] += -N3**2*detJ*r13*r33*wij + k += 1 + CAv[k] += -N3*N4*detJ*r13**2*wij + k += 1 + CAv[k] += -N3*N4*detJ*r13*r23*wij + k += 1 + CAv[k] += -N3*N4*detJ*r13*r33*wij + k += 1 + CAv[k] += -N1*N3*detJ*r13*r23*wij + k += 1 + CAv[k] += -N1*N3*detJ*r23**2*wij + k += 1 + CAv[k] += -N1*N3*detJ*r23*r33*wij + k += 1 + CAv[k] += -N2*N3*detJ*r13*r23*wij + k += 1 + CAv[k] += -N2*N3*detJ*r23**2*wij + k += 1 + CAv[k] += -N2*N3*detJ*r23*r33*wij + k += 1 + CAv[k] += -N3**2*detJ*r13*r23*wij + k += 1 + CAv[k] += -N3**2*detJ*r23**2*wij + k += 1 + CAv[k] += -N3**2*detJ*r23*r33*wij + k += 1 + CAv[k] += -N3*N4*detJ*r13*r23*wij + k += 1 + CAv[k] += -N3*N4*detJ*r23**2*wij + k += 1 + CAv[k] += -N3*N4*detJ*r23*r33*wij + k += 1 + CAv[k] += -N1*N3*detJ*r13*r33*wij + k += 1 + CAv[k] += -N1*N3*detJ*r23*r33*wij + k += 1 + CAv[k] += -N1*N3*detJ*r33**2*wij + k += 1 + CAv[k] += -N2*N3*detJ*r13*r33*wij + k += 1 + CAv[k] += -N2*N3*detJ*r23*r33*wij + k += 1 + CAv[k] += -N2*N3*detJ*r33**2*wij + k += 1 + CAv[k] += -N3**2*detJ*r13*r33*wij + k += 1 + CAv[k] += -N3**2*detJ*r23*r33*wij + k += 1 + CAv[k] += -N3**2*detJ*r33**2*wij + k += 1 + CAv[k] += -N3*N4*detJ*r13*r33*wij + k += 1 + CAv[k] += -N3*N4*detJ*r23*r33*wij + k += 1 + CAv[k] += -N3*N4*detJ*r33**2*wij + k += 1 + CAv[k] += -N1*N4*detJ*r13**2*wij + k += 1 + CAv[k] += -N1*N4*detJ*r13*r23*wij + k += 1 + CAv[k] += -N1*N4*detJ*r13*r33*wij + k += 1 + CAv[k] += -N2*N4*detJ*r13**2*wij + k += 1 + CAv[k] += -N2*N4*detJ*r13*r23*wij + k += 1 + CAv[k] += -N2*N4*detJ*r13*r33*wij + k += 1 + CAv[k] += -N3*N4*detJ*r13**2*wij + k += 1 + CAv[k] += -N3*N4*detJ*r13*r23*wij + k += 1 + CAv[k] += -N3*N4*detJ*r13*r33*wij + k += 1 + CAv[k] += -N4**2*detJ*r13**2*wij + k += 1 + CAv[k] += -N4**2*detJ*r13*r23*wij + k += 1 + CAv[k] += -N4**2*detJ*r13*r33*wij + k += 1 + CAv[k] += -N1*N4*detJ*r13*r23*wij + k += 1 + CAv[k] += -N1*N4*detJ*r23**2*wij + k += 1 + CAv[k] += -N1*N4*detJ*r23*r33*wij + k += 1 + CAv[k] += -N2*N4*detJ*r13*r23*wij + k += 1 + CAv[k] += -N2*N4*detJ*r23**2*wij + k += 1 + CAv[k] += -N2*N4*detJ*r23*r33*wij + k += 1 + CAv[k] += -N3*N4*detJ*r13*r23*wij + k += 1 + CAv[k] += -N3*N4*detJ*r23**2*wij + k += 1 + CAv[k] += -N3*N4*detJ*r23*r33*wij + k += 1 + CAv[k] += -N4**2*detJ*r13*r23*wij + k += 1 + CAv[k] += -N4**2*detJ*r23**2*wij + k += 1 + CAv[k] += -N4**2*detJ*r23*r33*wij + k += 1 + CAv[k] += -N1*N4*detJ*r13*r33*wij + k += 1 + CAv[k] += -N1*N4*detJ*r23*r33*wij + k += 1 + CAv[k] += -N1*N4*detJ*r33**2*wij + k += 1 + CAv[k] += -N2*N4*detJ*r13*r33*wij + k += 1 + CAv[k] += -N2*N4*detJ*r23*r33*wij + k += 1 + CAv[k] += -N2*N4*detJ*r33**2*wij + k += 1 + CAv[k] += -N3*N4*detJ*r13*r33*wij + k += 1 + CAv[k] += -N3*N4*detJ*r23*r33*wij + k += 1 + CAv[k] += -N3*N4*detJ*r33**2*wij + k += 1 + CAv[k] += -N4**2*detJ*r13*r33*wij + k += 1 + CAv[k] += -N4**2*detJ*r23*r33*wij + k += 1 + CAv[k] += -N4**2*detJ*r33**2*wij + From ab2aa398ca4b36d59d624b1581f6f5f99a053541 Mon Sep 17 00:00:00 2001 From: "Saullo G. P. Castro" Date: Wed, 11 Oct 2023 23:17:16 +0200 Subject: [PATCH 3/8] MAINT: updated scope and version --- setup.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 8ba222e..4850106 100644 --- a/setup.py +++ b/setup.py @@ -9,7 +9,7 @@ is_released = True -version = '0.3.24' +version = '0.4.0' def git_version(): @@ -192,6 +192,7 @@ def read(fname): 'static analysis', 'buckling', 'vibration', + 'panel flutter', 'structural dynamics', 'implicit time integration', 'explicit time integration', From ff90b13e9176df94c71d34db1875f036b7429bd1 Mon Sep 17 00:00:00 2001 From: "Saullo G. P. Castro" Date: Wed, 11 Oct 2023 23:19:02 +0200 Subject: [PATCH 4/8] MAINT: updated version number --- CITATION.cff | 6 +++--- README.md | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/CITATION.cff b/CITATION.cff index b260af7..e40bec7 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -16,7 +16,7 @@ abstract: >- Colab environment. message: >- Please, cite this software using: - Saullo G. P. Castro. (2023). General-purpose finite element solver based on Python and Cython (Version 0.3.24). Zenodo. DOI: https://doi.org/10.5281/zenodo.6573489. + Saullo G. P. Castro. (2023). General-purpose finite element solver based on Python and Cython (Version 0.4.0). Zenodo. DOI: https://doi.org/10.5281/zenodo.6573489. authors: - given-names: Saullo G. P. family-names: Castro @@ -30,5 +30,5 @@ identifiers: repository-code: 'https://github.com/saullocastro/pyfe3d' url: 'https://saullocastro.github.io/pyfe3d/' license: BSD-3-Clause -version: 0.3.24 -date-released: '2023-03-26' +version: 0.4.0 +date-released: '2023-10-11' diff --git a/README.md b/README.md index a72a7dc..645329f 100644 --- a/README.md +++ b/README.md @@ -22,7 +22,7 @@ in any platform, including the Google Colab environment. Citing this library ------------------- -Saullo G. P. Castro. (2023). General-purpose finite element solver based on Python and Cython (Version 0.3.24). Zenodo. DOI: https://doi.org/10.5281/zenodo.6573489. +Saullo G. P. Castro. (2023). General-purpose finite element solver based on Python and Cython (Version 0.4.0). Zenodo. DOI: https://doi.org/10.5281/zenodo.6573489. Documentation From 6232ee1cdc0ff3714d55b75b91eb3cab26ba4a0e Mon Sep 17 00:00:00 2001 From: "Saullo G. P. Castro" Date: Wed, 11 Oct 2023 23:19:35 +0200 Subject: [PATCH 5/8] TST: new test for piston theory --- tests/test_quad4r_piston_theory.py | 218 +++++++++++++++++++++++++++++ 1 file changed, 218 insertions(+) create mode 100644 tests/test_quad4r_piston_theory.py diff --git a/tests/test_quad4r_piston_theory.py b/tests/test_quad4r_piston_theory.py new file mode 100644 index 0000000..f0ff846 --- /dev/null +++ b/tests/test_quad4r_piston_theory.py @@ -0,0 +1,218 @@ +import sys +sys.path.append('..') + +import numpy as np +from numpy import isclose +from scipy.sparse.linalg import eigsh, eigs +from scipy.sparse import coo_matrix + +from pyfe3d.shellprop_utils import isotropic_plate +from pyfe3d import Quad4R, Quad4RData, Quad4RProbe, INT, DOUBLE, DOF + + +def test_quad4r_piston_theory(plot=False, refinement=1): + data = Quad4RData() + probe = Quad4RProbe() + nx = refinement*21 + ny = refinement*21 + + a = 0.8 + b = 0.5 + + E = 70.e9 # Pa + nu = 0.3 + + rho = 7.8e3 # kg/m3 + h = 0.0035 # m + + xtmp = np.linspace(0, a, nx) + ytmp = np.linspace(0, b, ny) + + dx = xtmp[1] - xtmp[0] + dy = ytmp[1] - ytmp[0] + + xmesh, ymesh = np.meshgrid(xtmp, ytmp) + ncoords = np.vstack((xmesh.T.flatten(), ymesh.T.flatten(), np.zeros_like(ymesh.T.flatten()))).T + + x = ncoords[:, 0] + y = ncoords[:, 1] + z = ncoords[:, 2] + ncoords_flatten = ncoords.flatten() + + inner = np.logical_not(isclose(x, 0) | isclose(x, a) | isclose(y, 0) | isclose(y, b)) + np.random.seed(20) + rdm = (-1 + 2*np.random.rand(x[inner].shape[0])) + np.random.seed(20) + rdm = (-1 + 2*np.random.rand(y[inner].shape[0])) + x[inner] += dx*rdm*0.4 + y[inner] += dy*rdm*0.4 + + nids = 1 + np.arange(ncoords.shape[0]) + nid_pos = dict(zip(nids, np.arange(len(nids)))) + nids_mesh = nids.reshape(nx, ny) + n1s = nids_mesh[:-1, :-1].flatten() + n2s = nids_mesh[1:, :-1].flatten() + n3s = nids_mesh[1:, 1:].flatten() + n4s = nids_mesh[:-1, 1:].flatten() + + num_elements = len(n1s) + print('num_elements', num_elements) + + KC0r = np.zeros(data.KC0_SPARSE_SIZE*num_elements, dtype=INT) + KC0c = np.zeros(data.KC0_SPARSE_SIZE*num_elements, dtype=INT) + KC0v = np.zeros(data.KC0_SPARSE_SIZE*num_elements, dtype=DOUBLE) + Mr = np.zeros(data.M_SPARSE_SIZE*num_elements, dtype=INT) + Mc = np.zeros(data.M_SPARSE_SIZE*num_elements, dtype=INT) + Mv = np.zeros(data.M_SPARSE_SIZE*num_elements, dtype=DOUBLE) + KA_betar = np.zeros(data.KA_BETA_SPARSE_SIZE*num_elements, dtype=INT) + KA_betac = np.zeros(data.KA_BETA_SPARSE_SIZE*num_elements, dtype=INT) + KA_betav = np.zeros(data.KA_BETA_SPARSE_SIZE*num_elements, dtype=DOUBLE) + KA_gammar = np.zeros(data.KA_GAMMA_SPARSE_SIZE*num_elements, dtype=INT) + KA_gammac = np.zeros(data.KA_GAMMA_SPARSE_SIZE*num_elements, dtype=INT) + KA_gammav = np.zeros(data.KA_GAMMA_SPARSE_SIZE*num_elements, dtype=DOUBLE) + CAr = np.zeros(data.CA_SPARSE_SIZE*num_elements, dtype=INT) + CAc = np.zeros(data.CA_SPARSE_SIZE*num_elements, dtype=INT) + CAv = np.zeros(data.CA_SPARSE_SIZE*num_elements, dtype=DOUBLE) + + N = DOF*nx*ny + + prop = isotropic_plate(thickness=h, E=E, nu=nu, calc_scf=True, rho=rho) + + quads = [] + init_k_KC0 = 0 + init_k_M = 0 + init_k_KA_beta = 0 + init_k_KA_gamma = 0 + init_k_CA = 0 + for n1, n2, n3, n4 in zip(n1s, n2s, n3s, n4s): + pos1 = nid_pos[n1] + pos2 = nid_pos[n2] + pos3 = nid_pos[n3] + pos4 = nid_pos[n4] + r1 = ncoords[pos1] + r2 = ncoords[pos2] + r3 = ncoords[pos3] + normal = np.cross(r2 - r1, r3 - r2)[2] + assert normal > 0 + quad = Quad4R(probe) + quad.n1 = n1 + quad.n2 = n2 + quad.n3 = n3 + quad.n4 = n4 + quad.c1 = DOF*nid_pos[n1] + quad.c2 = DOF*nid_pos[n2] + quad.c3 = DOF*nid_pos[n3] + quad.c4 = DOF*nid_pos[n4] + quad.init_k_KC0 = init_k_KC0 + quad.init_k_M = init_k_M + quad.init_k_KA_beta = init_k_KA_beta + quad.init_k_KA_gamma = init_k_KA_gamma + quad.init_k_CA = init_k_CA + quad.update_rotation_matrix(ncoords_flatten) + quad.update_probe_xe(ncoords_flatten) + quad.update_KC0(KC0r, KC0c, KC0v, prop) + quad.update_M(Mr, Mc, Mv, prop) + quad.update_KA_beta(KA_betar, KA_betac, KA_betav) + quad.update_KA_gamma(KA_gammar, KA_gammac, KA_gammav) + quad.update_CA(CAr, CAc, CAv) + quads.append(quad) + init_k_KC0 += data.KC0_SPARSE_SIZE + init_k_M += data.M_SPARSE_SIZE + init_k_KA_beta += data.KA_BETA_SPARSE_SIZE + init_k_KA_gamma += data.KA_GAMMA_SPARSE_SIZE + init_k_CA += data.CA_SPARSE_SIZE + + print('elements created') + + KC0 = coo_matrix((KC0v, (KC0r, KC0c)), shape=(N, N)).tocsc() + M = coo_matrix((Mv, (Mr, Mc)), shape=(N, N)).tocsc() + KA_beta = coo_matrix((KA_betav, (KA_betar, KA_betac)), shape=(N, N)).tocsc() + KA_gamma = coo_matrix((KA_gammav, (KA_gammar, KA_gammac)), shape=(N, N)).tocsc() + CA = coo_matrix((CAv, (CAr, CAc)), shape=(N, N)).tocsc() + + print('sparse KC0 and M created') + + bk = np.zeros(N, dtype=bool) + edges = np.isclose(x, 0.) | np.isclose(x, a) | np.isclose(y, 0) | np.isclose(y, b) + bk[0::DOF][edges] = True + bk[1::DOF][edges] = True + bk[2::DOF][edges] = True + + bu = ~bk + + Kuu = KC0[bu, :][:, bu] + Muu = M[bu, :][:, bu] + KA_betauu = KA_beta[bu, :][:, bu] + KA_gammauu = KA_gamma[bu, :][:, bu] + CAuu = CA[bu, :][:, bu] + + num_eigenvalues = 7 + print('eig solver begins') + # solves Ax = lambda M x + # we have Ax - lambda M x = 0, with lambda = omegan**2 + eigvals, eigvecsu = eigsh(A=Kuu, M=Muu, sigma=-1., which='LM', + k=num_eigenvalues, tol=1e-3) + print('eig solver end') + eigvecs = np.zeros((N, eigvecsu.shape[1]), dtype=float) + eigvecs[bu, :] = eigvecsu + omegan = eigvals**0.5 + + # panel flutter analysis + def MAC(mode1, mode2): + return (mode1@mode2)**2/((mode1@mode1)*(mode2@mode2)) + + MACmatrix = np.zeros((num_eigenvalues, num_eigenvalues)) + rho_air = 1.225 # kg/m^3 + v_sound = 343 # m/s + v_air = np.linspace(1.1*v_sound, 5*v_sound, 20) + Mach = v_air/v_sound + betas = rho_air*v_air**2/np.sqrt(Mach**2 - 1) + radius = 1. + gammas = betas/(2*radius*np.sqrt(Mach**2 - 1)) + mus = betas/(Mach*v_sound)*((Mach**2 - 2)/(Mach**2 - 1)) + mus*=0 # TODO quadractic eigenvalue problem to solve including damping + + omegan_vec = [] + for i, beta in enumerate(betas): + print('analysis i', i) + # solving generalized eigenvalue problem + Keffective = Kuu + beta*KA_betauu + gammas[i]*KA_gammauu + eigvals, eigvecsu = eigs(A=Keffective, M=Muu, + k=num_eigenvalues, which='LM', sigma=-1., tol=1e-6, v0=eigvecsu[:, 0]) + # NOTE with v0=eigvecsu[:, 0]) avoids fluctuations in adjacent + # solutions, while making it a bit slower + eigvecs = np.zeros((N, num_eigenvalues), dtype=float) + eigvecs[bu, :] = eigvecsu + + if i == 0: + eigvecs_ref = eigvecs + + corresp = [] + for j in range(num_eigenvalues): + for k in range(num_eigenvalues): + MACmatrix[j, k] = MAC(eigvecs_ref[:, j], eigvecs[:, k]) + if np.isclose(np.max(MACmatrix[j, :]), 1.): + corresp.append(np.argmax(MACmatrix[j, :])) + else: + corresp.append(j) + omegan_vec.append(eigvals[corresp]**0.5) + print(np.round(MACmatrix, 2)) + + eigvecs_ref = eigvecs[:, corresp].copy() + + + omegan_vec = np.array(omegan_vec) + + if plot: + import matplotlib + matplotlib.use('TkAgg') + import matplotlib.pyplot as plt + for i in range(num_eigenvalues): + plt.plot(Mach, omegan_vec[:, i]) + plt.ylabel('$\omega_n\ [rad/s]$') + plt.xlabel('Mach') + plt.show() + + +if __name__ == '__main__': + test_quad4r_piston_theory(plot=True, refinement=1) From 82aa197634b1764c1f11a2a16d3867b744227380 Mon Sep 17 00:00:00 2001 From: "Saullo G. P. Castro" Date: Wed, 11 Oct 2023 23:24:00 +0200 Subject: [PATCH 6/8] BLD: trying to fix new build --- pyproject.toml | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 pyproject.toml diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..7c5431c --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,2 @@ +[build-system] +requires = ["setuptools", "wheel", "numpy", "scipy", "alg3dpy", "cython"] From 1c9e39977fdf2f11f9d97f6a5928a663010a7e2e Mon Sep 17 00:00:00 2001 From: "Saullo G. P. Castro" Date: Thu, 12 Oct 2023 09:03:27 +0200 Subject: [PATCH 7/8] BUG: fixed casting of arrays in some systems + ENH: improved performance --- pyfe3d/shellprop.pyx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyfe3d/shellprop.pyx b/pyfe3d/shellprop.pyx index b3938e3..4311a79 100644 --- a/pyfe3d/shellprop.pyx +++ b/pyfe3d/shellprop.pyx @@ -1007,6 +1007,6 @@ cdef class GradABDE: self.gradEij[2, 0] = (mat.u6 + (-1)*mat.u7*lp.xiE1 + 0*lp.xiE2) # d(E11 E12 E16 E22 E26 E66) / d(xiE1, xiE2) - self.gradEij[0, 1:] = h*np.array([mat.u7, 0]) - self.gradEij[1, 1:] = h*np.array([0, -mat.u7]) - self.gradEij[2, 1:] = h*np.array([-mat.u7, 0]) + self.gradEij[0, 1] = h*mat.u7 + self.gradEij[1, 2] = h*(-mat.u7) + self.gradEij[2, 1] = h*(-mat.u7) From 738dd1d82238efce4d85bf4e1c416e9d404faa20 Mon Sep 17 00:00:00 2001 From: "Saullo G. P. Castro" Date: Thu, 12 Oct 2023 09:03:58 +0200 Subject: [PATCH 8/8] TST: more stable solver --- tests/test_tria3r_static_point_load.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/tests/test_tria3r_static_point_load.py b/tests/test_tria3r_static_point_load.py index 9ed17b8..e7635db 100644 --- a/tests/test_tria3r_static_point_load.py +++ b/tests/test_tria3r_static_point_load.py @@ -2,7 +2,7 @@ sys.path.append('..') import numpy as np -from scipy.sparse.linalg import cg +from scipy.sparse.linalg import spsolve from scipy.sparse import coo_matrix from pyfe3d.shellprop_utils import isotropic_plate @@ -14,6 +14,7 @@ def test_tria3r_static_point_load(plot=False, refinement=1): probe = Tria3RProbe() nx = 7*refinement ny = 11*refinement + if (nx % 2) == 0: nx += 1 if (ny % 2) == 0: @@ -110,25 +111,24 @@ def test_tria3r_static_point_load(plot=False, refinement=1): # applying boundary conditions # simply supported bk = np.zeros(N, dtype=bool) - check = np.isclose(x, 0.) | np.isclose(x, a) | np.isclose(y, 0) | np.isclose(y, b) - bk[2::DOF] = check - bk[0::DOF] = check - bk[1::DOF] = check + edges = np.isclose(x, 0.) | np.isclose(x, a) | np.isclose(y, 0) | np.isclose(y, b) + bk[2::DOF][edges] = True + bk[0::DOF][edges] = True + bk[1::DOF][edges] = True bu = ~bk # point load at center node f = np.zeros(N) fmid = 1. - check = np.isclose(x, a/2) & np.isclose(y, b/2) - f[2::DOF][check] = fmid + middle = np.isclose(x, a/2) & np.isclose(y, b/2) + f[2::DOF][middle] = fmid KC0uu = KC0[bu, :][:, bu] fu = f[bu] assert fu.sum() == fmid - uu, info = cg(KC0uu, fu, atol=1e-9) - assert info == 0 + uu = spsolve(KC0uu, fu) u = np.zeros(N) u[bu] = uu