C1 continuity functions available for quad and interval mesh_dolfinx?

dolfinx
1

Hello,

I want to know the available options for C1 continuous functions for quad and interval mesh elements. I looked for below 3 options, but not finding the implementation process.

  1. At present, the Hermite element is not supported for quad and interval mesh in dolfinx. Can you provide an example in github where you created the hermite element for triangle mesh, so that I can try to build the element manually for quad mesh (or if it is already available).

  2. The current Serendipity support is available but, it is not fully C1 continuous. The other option is using Isogeometric analysis (IGA), using NURBS, which I am not much familair for C1 contionuous elements. However, I found tIGAr which can be used for C1 continuity for quad/interval mesh. The elements are defined in 3D space which can be curved. Can you provide an example where this method is feasible?

  3. I want to deifne my functionspace of Hermite element with shape (3,) as
    dolfinx.fem.functionspace(mesh_l, basix.ufl.element( "HER", mesh_l.topology.cell_name(), deg, shape=(3, ))) how can i use mixed element of CG to obtain same C1 continuous behavior?

Thanks.

2

Hello FeniCS support team,

I am not able to work with C1 continuity functionspace in dolfinx. Can you suggest me a way to approach. I am completely struck?

Does anyone already worked for Hermite functions in interval mesh (and quad)?

3

These are just not supported at the moment in the library. If you provide more details on what kind of equation you want to solve we might provide some further help. There are workarounds to using C1 elements such as DG methods, relaxation using mixed elements…

1 Like
4

@Jeremy_Bleyer

Thanks for your response. I am providing a MWE where the unknown function w (displacement) and dw/dx (rotation) need to be continuous. The mesh generated is a general circle in 3D space.

The major part of issue is the var form where w and grad(w) need to be continuous both.

Edit: I find there is error in generating the example mesh. I apologize for that. I could create another example interval mesh.


import dolfinx
import basix
from dolfinx.fem import form, petsc, Function, functionspace, locate_dofs_topological, apply_lifting, set_bc
import numpy as np
from mpi4py import MPI
from ufl import Jacobian, as_vector, dot, SpatialCoordinate, as_tensor, Measure, Mesh
from ufl import TrialFunction, TestFunction,lhs, rhs, dx, dot, eq, grad
from dolfinx.io import gmshio
import ufl
theta = np.linspace(0, 2 * np.pi, 20)
y = 5 * np.cos(theta)
z = 5 * np.sin(theta)
points = np.array([[0, yi, zi] for yi, zi in zip(y, z)])
cell=[]
for k in range(len(points)):
        cell.append([k,k+1])   
elem = np.array(cell)
el= basix.ufl.element("S", "interval", 1, shape=(3, ))
mesh_l = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, elem, points, ufl.Mesh(el))# Define local frame
t = Jacobian(mesh_l)
t1 = as_vector([t[0, 0], t[1, 0], t[2, 0]])
e2=  t1/ sqrt(dot(t1, t1))
e1=  as_vector([1,0,0]) # Right Lay up
e3= cross(e1,e2)
e3=  e3/ sqrt(dot(e3, e3))   

V_l = dolfinx.fem.functionspace(mesh_l, basix.ufl.element(
    "S", mesh_l.topology.cell_name(), 2, shape=(3, )))

e1l,e2l,e3l=Function(V_l), Function(V_l), Function(V_l)

fexpr1=dolfinx.fem.Expression(e1,V_l.element.interpolation_points(), comm=MPI.COMM_WORLD)
e1l.interpolate(fexpr1) 

fexpr2=dolfinx.fem.Expression(e2,V_l.element.interpolation_points(), comm=MPI.COMM_WORLD)
e2l.interpolate(fexpr2) 

fexpr3=dolfinx.fem.Expression(e3,V_l.element.interpolation_points(), comm=MPI.COMM_WORLD)
e3l.interpolate(fexpr3) 

x, dx=SpatialCoordinate(mesh_l), Measure('dx')(domain=mesh_l)
e=[e1l,e2l,e3l]

k22= deri(e)
x11,O=local_grad(e[0],x[0]),local_grad(e[0],x[1])
x22,x32= local_grad(e[1],x[1]), local_grad(e[1],x[2])
y2,y3=  local_grad(e[2],x[1]), local_grad(e[2],x[2])      # rotation matrix terms
Eps2= as_tensor([(x11, O,  x[2], -x[1]),
                        (O, O,  O, O),
                        (O,x[1]*x32-x[2]*x22 , O, O),
                        (O,O,y3,y2),
                        (O,O,O,O), 
                        (O,-2-0.5*k22*(x[1]*x32-x[2]*x22),O,O)])
d22= as_ vector([x11, x22,x32)])
def G1(w):
        w_d2= [local_grad(e[1],w[i]) for i in range(3)]
        w_d22=[local_grad(e[1],w_d2[i]) for i in range(3)]
        w_22=[local_grad(d22,w[ii]) for ii in range(3)]

        G2=local_grad(e[1],w[1])+local_grad(e[1],w[2])
        G3=local_grad(e[1],w[0])
        G5=-ddot(w_22,d3)
        G6=-(k22)*0.5*G3
        return as_tensor([0,G2,G3,0,G5,G6])
    
def ksp_solve(A,F,V):
    w = Function(V)
    ksp = petsc4py.PETSc.KSP()
    ksp.create(comm=MPI.COMM_WORLD)
    ksp.setOperators(A)
    ksp.setType("preonly")
    ksp.getPC().setType("lu")
    ksp.getPC().setFactorSolverType("mumps")
    ksp.getPC().setFactorSetUpSolverType()
    ksp.getPC().getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)  # detect null pivots
    ksp.getPC().getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)  # do not compute null space again
    ksp.setFromOptions()
    ksp.solve(F, w.vector)
    w.vector.ghostUpdate(
        addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
    ksp.destroy()
    return w 
    
# Nullspace:
index_map = V_l.dofmap.index_map
nullspace_basis = [dolfinx.la.create_petsc_vector(index_map, V_l.dofmap.index_map_bs) for i in range(3)]
with ExitStack() as stack:
    vec_local = [stack.enter_context(xx.localForm()) for xx in nullspace_basis]
    basis = [np.asarray(xx) for xx in vec_local]
dofs = [V_l.sub(i).dofmap.list for i in range(3)]      # Dof indices for each subspace (x, y and z dofs)
for i in range(3):                                   # Build translational null space basis
    basis[i][dofs[i]] = 1.0
dolfinx.la.orthonormalize(nullspace_basis)           # Create vector space basis and orthogonalize
nullspace= petsc4py.PETSc.NullSpace().create(nullspace_basis, comm=MPI.COMM_WORLD)

# Solve Var form
dv, v_= TrialFunction(V_l), TestFunction(V_l)
F2 = sum([dot(dot(Stiff_mat,G1(dv)+Eps2[:,0]), G1(v_))*dx])   
a,b=form(lhs(F2)),form(rhs(F2))
A,F = assemble_matrix(a),petsc.assemble_vector(b)
F.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)
A.setNullSpace(nullspace) 
nullspace.remove(F)
w=ksp_solve(A,F,V_l) # required unknown function
# this requires grad(w)-> rotations to be continuous with continuity of w (displacements)

Thanks for the support. I was looking if I can amke something with mixed interpolatiopn for seperating w and grad(w) over different fuunctionspace . I don’t if that could be the way and how?

Thanks

5

Hi @Jeremy_Bleyer and FeniCS support team,

Greetings!

I have gone through biharmonic example for implementing C1 continuity using DG functionspace method of jump and averages. I have a general query related to my implementation,

If my bilinear form contains multiple terms of first derivatives (requiring C1 continuity) like
a(w, v) = \int_{\Omega} \sum_{i=1}^{N} \frac{d w}{d x} \cdot f_i(x_2, x_3) \, v \, dx. To implement continuity constraint of first derivatives, does the below updated bilinear form works,
\begin{align} a(w, v) = &\int_{\Omega} \sum_{i=1}^{N} \frac{d w}{d x} \cdot f_i(x_2, x_3) \, v \, dx \\ &- \int_{E} \left\langle \frac{d w}{d x} \right\rangle \, [\![ v ]\!] \, ds \\ &- \int_{E} [\![ w ]\!] \, \left\langle \frac{d v}{d x} \right\rangle \, ds \\ &+ \alpha \int_{E} \frac{1}{h_E} [\![ w ]\!] \, [\![ v ]\!] \, ds. \end{align}
Meaning, do we need to apply a pair of jump and average for first derivative only once, or option b, for all 5 first derivatives terms separately in my original weak form.

Any help or reference is greatly appreciated.

6

I’m not entirely sure what you’re asking. But perhaps this paper will help: View article. It offers a derivation of the DG formulation of the biharmonic problem. Further you can reduce this down to the C^0-IPG formulation by considering a C^0 conforming basis (e.g. Lagrange elements) which eliminates most of the terms, and should leave you with what’s in the dolfinx demos.

1 Like