Read in quadratic mesh to represent CG2 elements

mesh
21

Sure (sorry) here I have one:
import numpy as np
from dolfinx.fem import (
Function,
FunctionSpace,
locate_dofs_geometrical,
dirichletbc,
Constant,
assemble_scalar,
set_bc,
form,
)

import ufl
#Librairy for p-FE
import basix
import basix.ufl_wrapper
from mpi4py import MPI
comm = MPI.COMM_WORLD
Lx = 2
Nb_elem= 200
from dolfinx.mesh import create_interval
mesh = create_interval(comm, Nb_elem, [0,Lx])

element = basix.create_element(basix.ElementFamily.P, basix.cell.string_to_type(mesh.ufl_cell().cellname()), 4, basix.LagrangeVariant.legendre, True)
element_u = basix.ufl_wrapper.BasixElement(element)
V_u = FunctionSpace(mesh, element_u)

u = Function(V_u, name=“Displacement”)
Uimp_ux = Function(V_u, name=“Boundary Displacement”)
zero_ux = Function(V_u, name=“Null Boundary Displacement”)

Boundary definitions dofs

dofs_u_left = locate_dofs_geometrical(V_u, lambda x: np.isclose(x[0], 0.0))
dofs_u_right = locate_dofs_geometrical(V_u, lambda x: np.isclose(x[0], Lx))

Set Bcs Function

zero_ux.interpolate(lambda x: (np.zeros_like(x[0])))
Uimp_ux.interpolate(lambda x: (np.ones_like(x[0])))

bcs_u = [
dirichletbc(zero_ux, dofs_u_left),
dirichletbc(Uimp_ux, dofs_u_right),
]

Thanks

22

Please use 3x` encapsulation around code, i.e.

```python
import dolfinx
#add code here
```
23

Sure

from dolfinx.fem import (
Function,
FunctionSpace,
locate_dofs_geometrical,
dirichletbc,
Constant,
assemble_scalar,
set_bc,
form,
)

import ufl
#Librairy for p-FE
import basix
import basix.ufl_wrapper
from mpi4py import MPI
comm = MPI.COMM_WORLD
Lx = 2
Nb_elem= 200
from dolfinx.mesh import create_interval
mesh = create_interval(comm, Nb_elem, [0,Lx])

element = basix.create_element(basix.ElementFamily.P, basix.cell.string_to_type(mesh.ufl_cell().cellname()), 4, basix.LagrangeVariant.legendre, True)
element_u = basix.ufl_wrapper.BasixElement(element)
V_u = FunctionSpace(mesh, element_u)

u = Function(V_u, name=“Displacement”)
Uimp_ux = Function(V_u, name=“Boundary Displacement”)
zero_ux = Function(V_u, name=“Null Boundary Displacement”) 

dofs_u_left = locate_dofs_geometrical(V_u, lambda x: np.isclose(x[0], 0.0))
dofs_u_right = locate_dofs_geometrical(V_u, lambda x: np.isclose(x[0], Lx))

zero_ux.interpolate(lambda x: (np.zeros_like(x[0])))
Uimp_ux.interpolate(lambda x: (np.ones_like(x[0])))

bcs_u = [
dirichletbc(zero_ux, dofs_u_left),
dirichletbc(Uimp_ux, dofs_u_right),
] ```
24

You cannot use strong boundary conditions with Legendre elements, as:

  1. They are discontinuous, all dofs belong to the interior of the cell
  2. The interpolation matrix is not an identity operator (they are integral moments).

You would have to use DG/Nitsche type enforcement of bcs.

25

Hello, following your answer I would have a more general question. I am doing a PhD and I plan in the future to use CG Legendre high-order shape functions (until p=8) in multi-dimensional problems.

I would like to know if it is possible in FEniCSx to create these shape functions (although they are not available today according to my understanding) and also to manage the assembly of the matrixes that is complicated especially with tetrahedron.

Thank you for your help

26

Im not sure it is sensible to use CG with a Legendre basis, see: Legendre polynomial as basis for finite element method - Mathematics Stack Exchange
and
finite element - 1D Discontinuous Galerkin - Lagrange vs Legendre Basis - Computational Science Stack Exchange

Could you elaborate on why you want to use this?

27

Sure, so the final idea is to import a xdmf file of a 3D mesh of a long bone with tetrahedrons. Then we would like to solve an elasticity problem as a first step (so this should include the process of assembling the stiffness matrix).
Finally we will consider to solve a Phase Field problem (weakly coupled problem elasticity + damage) in order to study the bone fracture.

Thank you again

28

To me this still has to be a DG formulation, Since Legendre polynomials as a basis does not associate its functional with a vertex, Edge or facet.

29

Actually, I may not have explain well, I meant using integrated Legendre polynomials as shape functions (Szabo and Babuska shape functions) that are hierarchical. They demonstrated good convergence properties (in terms of DOFs compared to Lagrange shape functions and CPU time whether the hierarchical property is used)

30

I’ve seen these referred to as Lobatto polynomials before. These can be defined in FEniCSx using Basix’s custom element functionality. I’ve had some success implementing them for interval, quadrilateral, an hexahedron cells, but not look at these elements on triangles and tets yet.

31

Thank you very much for your answer.
I checked and it seems to be what I am looking for in the sense it is CG and also hierarchical.

Thanks again