On Sat, Mar 23, 2024 at 8:03 AM Gong Yujie <yc17...@connect.um.edu.mo> wrote:
> Dear PETSc group, I'm reading the DS part for the discretization start > from SNES ex17. c which is a demo for solving linear elasticity problem. I > have two questions for the details. The first question is for the residual > function. Is the residual > ZjQcmQRYFpfptBannerStart > This Message Is From an External Sender > This message came from outside your organization. > > ZjQcmQRYFpfptBannerEnd > Dear PETSc group, > > I'm reading the DS part for the discretization start from SNES ex17.c > which is a demo for solving linear elasticity problem. I have two questions > for the details. > > The first question is for the residual function. Is the residual > calculated as this? The dot product is a little weird because of the > dimension of the result. > Here \sigma is the stress tensor, \phi_i is the test function for the i-th > function (Linear elasticity in 3D contains three equations). > The stress term in the momentum equation is (-div sigma) . psi = d_i sigma_{ij} psi_j which is then integrated by parts sigma_{ij} d_i psi_j This is linear isotropic elasticity, so sigma = \mu (d_i u_j + d_j u_i) + \lambda \delta_{ij} sigma_{kk} In PETSc, we phrase the term in the weak form as f^1_{ij} d_i psi_j so f1[i * dim + j] below is sigma_{ij}, from line 298 of ex17.c for (PetscInt c = 0; c < Nc; ++c) { for (PetscInt d = 0; d < dim; ++d) { f1[c * dim + d] += mu * (u_x[c * dim + d] + u_x[d * dim + c]); f1[c * dim + c] += lambda * u_x[d * dim + d]; } } > The second question is how to derive the Jacobian of the system (line 330 > in ex17.c). As shown in the PetscDSSetJacobian, we need to provide function > g3() which I think is a 4th-order tensor with size 3*3*3*3 in this linear > elasticity case. I'm not sure how to get it. Are there any references on > how to get this Jacobian? > The Jacobian indices are shown here: https://urldefense.us/v3/__https://petsc.org/main/manualpages/FE/PetscFEIntegrateJacobian/__;!!G_uCfscf7eWS!eNdXZ9QpzXPTMCtLneeS7td5OYVeKmhdOrw_fswDze2u7v2JaO7kBmFTVakyHDmWJOCHjf-0GrqGc53QKQec$ where the g3 term is \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q) To get the Jacobian, we use u = \sum_i u_i psi_i, where psi_i is a vector, and then differentiate the expression with respect to the coefficient u_i. Since the operator is already linear,this is just matching indices for (PetscInt c = 0; c < Nc; ++c) { for (PetscInt d = 0; d < dim; ++d) { g3[((c * Nc + c) * dim + d) * dim + d] += mu; g3[((c * Nc + d) * dim + d) * dim + c] += mu; g3[((c * Nc + d) * dim + c) * dim + d] += lambda; } } Take the first mu term, mu (d_c u_d ) (d_c psi_d). We know that fc == gc and df == dg, so we get g3[((c * Nc + c) * dim + d) * dim + d] += mu; For the second term, we transpose grad u, so fc == dg and gc == df, g3[((c * Nc + d) * dim + d) * dim + c] += mu; Finally, for the lambda term, fc == df and gc == dg because we are matching terms in the sum g3[((c * Nc + d) * dim + c) * dim + d] += lambda; Thanks, Matt > I've checked about the comment before this Jacobian function (line 330 in > ex17.c) but don't know how to get this. > > Thanks in advance! > > Best Regards, > Yujie > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!eNdXZ9QpzXPTMCtLneeS7td5OYVeKmhdOrw_fswDze2u7v2JaO7kBmFTVakyHDmWJOCHjf-0GrqGc4DkrV1c$ <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!eNdXZ9QpzXPTMCtLneeS7td5OYVeKmhdOrw_fswDze2u7v2JaO7kBmFTVakyHDmWJOCHjf-0GrqGcwApR5ek$ >