Stiffness matrix

[liuyuanhao]

Dear canitesc,
I want to learn your code now.But I meet a little difficult.May you tell me that you use which rules to build the stiffness matrix in the CantileverBeam.m?Do you have the paper or math formula about the implemention of the stiffness matrix?
Thank you so much!
Deeply grateful,
Jim

[canitesc]

Hi,
I think the assembly procedure is fairly standard, although it is slightly more complicated because of the IGA. The formulas and the original MATLAB FEM code on which this is based can be found e.g. here https://www.math.purdue.edu/~caiz/math615/matlab_fem.pdf The main difference is that the spline basis functions are computed using the Bezier extraction operator, some description of the method is given in this paper: https://doi.org/10.1002/nme.2968

[liuyuanhao]

Dear canitesc,
Thanks for your previous help,I want to ask you a question.Can this code package contain functions to compute the second derivative of PHT splines?
Deeply grateful,
Jim

[canitesc]

Hi Jim,
As Somdatta mentioned in somdattagoswami/IGAPack-PhaseField#2 , the code for second derivative are included in her repository. Actually, since the PHT-splines are based on the Bernstein-Bézier form, you mainly just need the second derivatives of the Bernstein polynomials.
-Cosmin

[liuyuanhao]

Dear canitesc,
Thanks for your previous help,I want to ask you a question.
In the code --assembleGalerkinSysGIFTMP.m,why the dimension of the dRdx is different from the dimension of the gauss_coord_x?Can I find the coordinate of the control points of the derivatives such as dRdx,dRdy?
Deeply grateful,
Jim

[canitesc]

The dRdx and dRdy arrays contain the derivatives of all the basis functions for a given Gauss point and element. The derivatives of the Bernstein polynomials at the Gauss points on the reference element [-1,1]^d, where d is the dimension of the space, are precomputed in dBdu and dBdv at the top of the function since we need to use them many times in each element. Then, for d=2, the derivatives of the basis functions are computed as

dRdx = (PHTelem{patchIndex}(i).C)*squeeze(dBdu(ii,jj,:));
dRdy = (PHTelem{patchIndex}(i).C)*squeeze(dBdv(ii,jj,:));

where i is the element index, ii and jj are the Gauss point indices in the u and v direction per element and (PHTelem{patchIndex}(i).C) is a matrix with m rows and (p+1)*(q+1) columns where each row corresponds to a basis function and stores its Bézier coefficients. To compute the second derivatives of the basis functions, you need to replace dBdu and dBdv with arrays that store the second derivatives of the Bernstein polynomials.

[liuyuanhao]

Thank you so much. And may ask a question about the use of parameter p and q,what are they used for？And which function is used for generating control points？ Deeply Grateful, Jim

…

---Original--- From: "Cosmin Anitescu"<notifications@github.com&gt; Date: Wed, Feb 24, 2021 19:01 PM To: "canitesc/IGAPack"<IGAPack@noreply.github.com&gt;; Cc: "liuyuanhao"<77428159@qq.com&gt;;"Author"<author@noreply.github.com&gt;; Subject: Re: [canitesc/IGAPack] Stiffness matrix (#4) The dRdx and dRdy arrays contain the derivatives of all the basis functions for a given Gauss point and element. The derivatives of the Bernstein polynomials at the Gauss points on the reference element [-1,1]^d, where d is the dimension of the space, are precomputed in dBdu and dBdv at the top of the function since we need to use them many times in each element. Then, for d=2, the derivatives of the basis functions are computed as dRdx = (PHTelem{patchIndex}(i).C)*squeeze(dBdu(ii,jj,:)); dRdy = (PHTelem{patchIndex}(i).C)*squeeze(dBdv(ii,jj,:)); where i is the element index, ii and jj are the Gauss point indices in the u and v direction per element and (PHTelem{patchIndex}(i).C) is a matrix with m rows and (p+1)*(q+1) columns where each row corresponds to a basis function and stores its Bézier coefficients. To compute the second derivatives of the basis functions, you need to replace dBdu and dBdv with arrays that store the second derivatives of the Bernstein polynomials. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or unsubscribe.

[canitesc]

The parameters p and q are the polynomial degrees of the approximation in the u and v direction, although right now I think the code only works when p=q>=3. The control points for the geometry in the GIFT method are defined the subroutine init2DGeometryGIFTMP. Note that in the GIFT method, the control points are fixed as the mesh is refined. In the isoparametric (or isogeometric) method, for the examples where the filename ends in "iso", the control points are recalculated by the refinement subroutine and the polynomial degree in this code is fixed at cubics ( p=q=3).

[liuyuanhao]

Thank you so much.I totally understand it.
And may ask a question about the thin shell example?If I want to implement the thin shell(kirchhoff-love shell) such as cantilever beam and Scordelis-Lo Shell problem,(with almost no thickness)can I use the pack to do this?Which parameters should I change?
Cannot thank you enough.
Deeply Grateful,
Jim