Issue in integration of stifness matrix for basis degree and order-quad

[sajadj11]

Hi,

I'm working with Nutils and previously had a composite plate with basis degree 5 and 18 elements, which worked for calculating stresses. Now, I need to introduce a circular hole using topo = topo0.trim(np.linalg.norm(shifted_geom) - radius, maxrefine=maxrefine). However, I'm encountering an issue when calculating the stiffness matrix. With basis degree p=5 and quadrature order 2p+1, I get the warning:

UserWarning: inexact integration for polynomial of degree 11 warnings.warn('inexact integration for polynomial of degree {}'.format(degree)) This warning disappears only when I reduce p to 2 or 3 and set the quadrature order to 2p. I need to maintain p=5 for my calculations. Can you please help me resolve this issue?

Thank you!

[gertjanvanzwieten]

Trimming creates an internal triangulation for quadrature, and on triangle elements we only implemented quadrature schemes up to order 6. Hence the warning: any polynomial of degree > 6 will not be integrated exactly. However the quadrature error is probably small, considering that the integrand is likely to be smooth over the tiny triangles that are formed against the trimming boundary. The full degree still applies to the untrimmed, square elements, as well as the non-triangulated parts of the element subdivisions.

Higher order schemes for triangles certainly exist, for instance those published by Dunavant ("high degree efficient symmetrical Gaussian quadrature rules for the triangle", 1985) and it makes sense to implement those, so we'll keep this issue open until we get to it. But until then this is a warning that you can probably safely ignore.

[sajadj11]

Hi,

Thank you very much for your previous response. I have another issue that related to the previous problem I had.

I am working with the following part of my code to build the second derivative and first derivative matrices:

# Building second derivative matrix 
ns.N_njk = 'd(d(basis_n, X_j), X_k)'

# Building first derivative matrix
ns.Y_nj = 'd(basis_n, X_j)'

ns.R1 = np.stack([
    np.stack([ns.Y[0,0], ns.zero, ns.zero], axis=0),
    np.stack([ns.zero, ns.Y[0,1], ns.zero], axis=0),
    np.stack([ns.Y[0,1], ns.Y[0,0], ns.zero], axis=0),
    np.stack([ns.zero, ns.zero, ns.N[0,0,0]], axis=0),
    np.stack([ns.zero, ns.zero, ns.N[0,1,1]], axis=0),
    np.stack([ns.zero, ns.zero, 2*ns.N[0,0,1]], axis=0)
], axis=0)
R = ns.R1
for i in range(1, len(basis)):
    stacked_rows = np.stack([
        np.stack([ns.Y[i, 0], ns.zero, ns.zero]),
        np.stack([ns.zero, ns.Y[i, 1], ns.zero]),
        np.stack([ns.Y[i, 1], ns.Y[i, 0], ns.zero]),
        np.stack([ns.zero, ns.zero, ns.N[i, 0, 0]]),
        np.stack([ns.zero, ns.zero, ns.N[i, 1, 1]]),
        np.stack([ns.zero, ns.zero, 2 * ns.N[i, 0, 1]])
    ])
    R = np.concatenate((R, stacked_rows), axis=1)
ns.R = R
# Combined stiffness matrix
ns.M = np.concatenate([np.concatenate([ns.A, ns.B], axis=1), 
                       np.concatenate([ns.B, ns.D], axis=1)], axis=0)
K = topo.integral(ns.eval_pj('(R_kp M_kq R_qj) d:X'), degree=order_quad).eval()

When I increase the number of elements, I encounter a segmentation fault (core dumped) error in K. I assumed increasing the number of elements would improve accuracy for my problem, but this issue is blocking my progress. Can you please help me resolve this?

Best Regards

[joostvanzwieten]

I think the problem is with the definition of R. The length of the basis is proportional to the number of elements. By appending a row to R for each basis function, you're creating a data structure with a memory foot print that is proportional to the number of elements as well. Instead, you can use tensor operations to build R. Example for Nutils 8 and up:

# Building second derivative matrix
ns.N_njk = 'd(d(basis_n, X_j), X_k)'

# Building first derivative matrix
ns.Y_nj = 'd(basis_n, X_j)'

ns.δ6 = function.eye(6)
ns.δ3 = function.eye(3)
ns.R_inj = 'δ6_i0 δ3_j0 Y_n0 + δ6_i1 δ3_j1 Y_n1 + δ6_i2 δ3_j0 Y_n1 + δ6_i2 δ3_j1 Y_n0 + δ3_j2 (δ6_i3 N_n00 + δ6_i4 N_n11 + 2 δ6_i5 N_n01)'
ns.R = np.reshape(ns.R, (6, -1))

# Combined stiffness matrix
ns.M = np.concatenate([np.concatenate([ns.A, ns.B], axis=1),
                       np.concatenate([ns.B, ns.D], axis=1)], axis=0)
K = topo.integral(ns.eval_pj('(R_kp M_kq R_qj) d:X'), degree=order_quad).eval()

It looks like you're building some kind of jacobian by hand. Are you aware that Nutils can do this for you using function.derivative?

[sajadj11]

I think the problem is with the definition of R. The length of the basis is proportional to the number of elements. By appending a row to R for each basis function, you're creating a data structure with a memory foot print that is proportional to the number of elements as well. Instead, you can use tensor operations to build R. Example for Nutils 8 and up:

# Building second derivative matrix
ns.N_njk = 'd(d(basis_n, X_j), X_k)'

# Building first derivative matrix
ns.Y_nj = 'd(basis_n, X_j)'

ns.δ6 = function.eye(6)
ns.δ3 = function.eye(3)
ns.R_inj = 'δ6_i0 δ3_j0 Y_n0 + δ6_i1 δ3_j1 Y_n1 + δ6_i2 δ3_j0 Y_n1 + δ6_i2 δ3_j1 Y_n0 + δ3_j2 (δ6_i3 N_n00 + δ6_i4 N_n11 + 2 δ6_i5 N_n01)'
ns.R = np.reshape(ns.R, (6, -1))

# Combined stiffness matrix
ns.M = np.concatenate([np.concatenate([ns.A, ns.B], axis=1),
                       np.concatenate([ns.B, ns.D], axis=1)], axis=0)
K = topo.integral(ns.eval_pj('(R_kp M_kq R_qj) d:X'), degree=order_quad).eval()

It looks like you're building some kind of jacobian by hand. Are you aware that Nutils can do this for you using function.derivative?

Hi

Thank you so much for your previous response. Your guidance on redefining the matrix R and using tensor operations has resolved the segmentation fault error I was encountering.

I truly appreciate your expertise and support.