I am trying to change the elements in the mesh of the hyperelasticity demo to quadratic elements. What I would like to do is to change the cells in the mesh to square cells. The demo I’m referring to is here: Hyperelasticity — DOLFIN documentation
What I’ve tried to do is to change the mesh to:
mesh = UnitCubeMesh.create(24,16,16, CellType.Type.quadrilateral)
instead of:
mesh = UnitCubeMesh(24, 16, 16)
This is the error that I am getting:
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
It sounds like you mean “quadrilateral” rather than “quadratic” (where the latter would usually refer to the polynomial degree of the shape functions used, not the shape of the elements). The 3D counterpart to quadrilateral elements are hexahedron elements. You can create a structure hexahedron mesh with
As a sidenote to @kamensky’s answer, if you want to use quadrilateral or hexahedral elements, it is strongly adviced to use dolfinx, as the support for non-simplex elements in dolfin is quite limited. See for instance my example of hyperelasticity: Hyperelasticity — FEniCSx tutorial
Which uses hexahedral elements
Thank you very much, I was confused between quadratic and quadrilateral elements and this helped solve the problem. As an additional question, would it be necessary to use the option of ‘quadrature_degree’ for the demo to run well?
By default, UFL tries to estimate an appropriate number of quadrature points to integrate forms with, based on what operations appear in each integrand. However, in complicated problems like finite strain hyperelasticity, this automatic estimate is typically way too conservative (i.e., using too many points), and leads to unnecessary computational cost. Thus, it’s most efficient to manually set the quadrature degree, e.g.,
(and similarly for ds and/or dS, if applicable). The appropriate quadrature degree is problem-dependent, but I often use twice the polynomial degree of shape functions as a reasonable default, i.e., enough to integrate products of shape functions, as found in mass matrices. In static elastic problems you can go lower, since you’re mainly interested in products of derivatives of shape functions. For example, using linear elasticity as a guide, with linear simplicial elements, you have constant strain on each element, so a quadrature degree of zero (i.e., one point per element) is sufficient. For bilinear/trilinear quad/hex elements, though, there are some quadratic terms in the shape functions, so you’d want to use at least a quadrature degree of one (i.e., a tensor-product quadrature rule with two points in each direction), since one-point quadrature would leave low-energy “hourglass” modes in the system.