Weak C1 Continuity in 2D (quad) mesh using CDG

General
1

Dear FeniCS team,

Greetings!
I want to implement C1 continuity using C0 elements in 2D mesh (quad). In the Biharmonic equation — DOLFINx 0.7.2 documentation, C1 continuity is ensured along normal derivative using Penalty method (CDG).
My query is for maintaining C1 continuity, we need to constraint both normal and tangent derivative along the interface (edge), whereas only normal derivative was constrained.

For my case of anisotropic material, I also need C1 continuity among tangent, where I cannot use inbuilt avg and jump functions, and n(+) and n(-) feature.

To implement manually, I thought to generate tangent vector like n=FacetNormal(mesh), but (n(+) and n(-) feature will be missing still. Can you guide me how I can implement the weak tangent continuity and use inbuilt features?

Note: For 1D mesh, C1 continuity is working for me.
Thanks

2

You can just define your own jump and average operators. E.g. for the tangent jump

As used in a C0-IPG scheme I have in a preprint here (eq 86). Or in a simpler context of a Maxwell operator here Section 4.6.

1 Like
3

@nate Thanks for the response. My 2D mesh is defined in a 3D space. Using dg_cross, works only for planar domain. I realize, it is generating the tangent vector by taking cross-product and use it in cross_jump.

I can use normal to element (3rd vector), by cross-product of Jacobian(mesh) vectors. This normal can be used to make cross-product for determining facet tangent and perform what dg_cross is expected to do. I apologize for anything that I miss in interpretation.

4

Sure, you can do whatever you like. I’m just emphasising that you can implement these operators yourself to satisfy your needs.

5

@nate

The Biharmonic example ensures continuity of first derivative along normal to the interface. I want to make weak continuity of cross-derivatives.

The MWE:

def cross_grad(dv):
    return as_tensor([((dv[0].dx(0)).dx(1), (dv[1].dx(0)).dx(1),(dv[2].dx(0)).dx(1)),
                     ((dv[0].dx(1)).dx(2), (dv[1].dx(1)).dx(2),(dv[2].dx(1)).dx(2)), 
                     ((dv[0].dx(0)).dx(2), (dv[1].dx(0)).dx(2),(dv[2].dx(0)).dx(2))]).T

h = CellDiameter(mesh)
h_avg = (h('+') + h('-')) / 2.0
dS=Measure('dS')(domain=mesh)
n = FacetNormal(mesh)
V= dolfinx.fem.functionspace(mesh, basix.ufl.element(
    "S", mesh.topology.cell_name(), 2, shape=(3, )))
dv,v_=TrialFunction(V), TestFunction(V)
dx=Measure('dx')(domain=mesh)
# Weak continuity bilinear terms of cross-derivatives:
cross_C0-IPM= - inner(avg(div(cross_grad(dv))), jump(cross_grad(v_), n))*dS \
    - inner(jump(cross_grad(dv), n), avg(div(cross_grad(v_))))*dS \
    + (alpha/h_avg**2)*inner(jump(cross_grad(dv),n), jump(cross_grad(v_),n))*dS

The above cross_C0-IPM doesn’t provide weak continuity of derivatives. My mesh is defined in 3D space with dv.ufl_shape=3. I tried to form a 3X3 tensor similar to grad(dv) and grad(v_), bt, it doesn’t work. Can you please comment or give feedback for my approach, it would be greatly helpful.

Thanks a lot!

6

If you want both the normal and tangential components, can’t you just add the entire jump?

inner(jump(grad(v_)),jump(grad(dv)))*dS

Edit: sorry, I wrote before reading the full thread. Although I don’t quite understand why the above wouldn’t work for you.

I am confused about your question though. You’re interested in C1 continuity, but in your latest MWE, you talk about cross derivatives, and have second derivatives in your operator. C1 concerns continuity of the first derivatives, so where are those 2nd derivatives coming from?

What is the math behind the operator that you’re trying to construct?

7

Yes, you are right. I was going in wrong direction to compute double derivative continuity. I came
back to first derivative continuity. In

,for all components of first derivative (grad(u) tensor to be continuous along interface, tangent facet and FacetNormal(mesh)XFacetTangent(mesh) directions to be used in addition to FacetNormal(mesh). On using

it directly was a good way, but, that doesn’t work. Only FacetNormal(mesh) provide partial support. The mesh is cylindrical shell defined in 3D space, where derivative has 3 components along x,y and z directions. Thanks a lot for your feedback. Any help is greatly appreciated.