Traction Separation Law in Dolfinx

General
1

All,

I am trying to leverage the work from @bleyerj from Cohesive zone modelling of debonding and bulk fracture to approach traction separation laws using dolfinx. I have worked to replicate the original code in the updated framework but am facing numerous challenges in its implementation, primarily in selecting appropriate function spaces and solver settings.

Has anyone done similar work in dolfinx?

2

Hi, sorry I did not have time yet to port this demo to dolfinx. I will try to do it in the near future though…
Can you be more specific on what is your exact issue ?

1 Like
3

@bleyerj , Thanks for the reply. My intent was just to ping the community, happy to provide specific details if you are interested.

Selection of Function Spaces
I have done little work with quadrature elements and spaces other than CG/DG. As per your comment in the above example

The ideal solution would be here to use Quadrature elements on the mesh facets, which are unfortunately not available in FEniCS. We will instead use a Discontinuous Lagrange Trace (also known as HDiv Trace ) space which represents polynomials living on the facet of the meshes and which are discontinuous at the vertices.

The issue hence is deciding the element family for V_f. I see from hdiv-trace-elements that DGT is not supported yet in Dolfinx. I am not familiar with using Quadrature elements in Dolfinx though I am trying to learn. Any advice is appreciated.

Thanks,

-Sven

4

Hi. Is there any update on this topic (CZM in FEniCSx)? I’m also interested. Thank you.

5

A first alternative is to update the code by Jeremy Bleyer using CR (Cruzier-Raviart) element for the jump field, with somenthing like

    "CR", msh.basix_cell(), surface_degree, discontinuous=False
)

I have done somenthing few months ago using this.

We started working on a better version with @bleyerj and @dokken two weeks ago. We should end up with somenthing hopefully for the end of the year. The plan is to use facet meshes and formulate the problem as a convex optimization problem, using Mosek as a solver.

1 Like
6

Excellent, thank you @cmaurini. I will keep exploring this. Eventually I aim at coupling it with some contact mechanics & friction to model a tool / obstacle, in the spirit of a cutting tool.

7

How is it going, do you think you will finish before the end of this year?

I want to define a hertzian-model between two cylinders and I think this would be a great help. Thank you in advance for pushing this project forward!

8

Hi, see here for the DG formulation eveywhere
https://bleyerj.github.io/comet-fenicsx/tours/interfaces/intrinsic_czm/intrinsic_czm.html
and here for a mixed CG formulation + a CZM at the interface
https://bleyerj.github.io/comet-fenicsx/tours/interfaces/czm_interface_only/czm_interface_only.html

1 Like
9

For contact problems you can check: Coupling PDEs of multiple dimensions — FEniCS Workshop

10

@bleyerj Thank you for sharing your work with the community. Very excited to dive into this formulation.

1 Like
11

Thank you for the tutorials!
One question regarding the mixed CG formulation:
For a non linear problem rather than linear problem, could the setup be somehow done in the spirit of:
problem = fem.petsc.NonlinearProblem(Res_u, [u1, u2], bc, ufl.derivative(Res_u, [u1, u2], [du1, du2]))?
Otherwise, how should the problem be implemented?

12

You can use the blocked solver from scifem: scifem/tests/test_blocked_newton_solver.py at main · scientificcomputing/scifem · GitHub
or hopefully soon the blocked/nest version with PETSc: SNES solver interface for standard, block and nest systems by jorgensd · Pull Request #3641 · FEniCS/dolfinx · GitHub

13

Thank you.

I have two additional questions:

1) Regarding the application of Dirichlet BC in the same tutorial (mixed CG):

Dirichlet BC are applied (for example for the left dofs) with left_dofs = fem.locate_dofs_topological(V1, fdim, subdomain1_facet_tags.find(1)), where the last argument refers to previously defined tags, also as in Multiphysics: Solving PDEs on subdomains — FEniCS Workshop.
My question: if I work with a mesh imported without identifiers (Physical groups, etc) and I rather want to locate the facets with a function, e.g.,

def bnd_encastre(x1):
    return np.isclose(...)

How should I proceed?

  • In the previous function, should x1 be set as x1 = ufl.SpatialCoordinate(subdomain1), with subdomain1 the submesh?
  • Should I locate the entities with dolfinx.mesh.locate_entities and then pass them to dolfinx.fem.locate_dofs_topological?
  • If yes, should it be done on the submesh (subdomain1) or on the parent mesh (mesh)?

Just to illustrate what I mean:

bnd_encastre_facets0_sub1 = dolfinx.mesh.locate_entities(subdomain1, mesh.topology.dim - 1, bnd_encastre)
bnd_encastre_dofs0_sub1 = dolfinx.fem.locate_dofs_topological(V1, mesh.topology.dim - 1, bnd_encastre_facets0_sub1)

bc_boundar_sub1 = dolfinx.fem.dirichletbc(np.array((0,) * mesh.geometry.dim, dtype=default_scalar_type), bnd_encastre_dofs0_sub1, V1)

Note that V1 refers to a function space on the submesh, the same as it is done in the tutorial.

2) I want to define a Function and interpolate to it an expression that depends i) on: the displacement (on each of the subdomains, u1, u2) and spatial position (x1, x2).

  • Should I define two functions (for each of the subdomains, e.g., v1 and v2), respectively from the function spaces on the submerges V1 and V2, and then do something similar to:
v1_expr = dolfinx.fem.Expression(function(x1,u1), V1.element.interpolation_points())
v1.interpolate(v1_expr)
v2_expr = dolfinx.fem.Expression(function(x2,u2), V2.element.interpolation_points())
v2.interpolate(v2_expr)

with function(x1, u1) denotes any generic expression defined previously.

I am bumping into errors in an implementation and I want to double check and understand the concepts before.

Thank you in advance.

14

It should be done of the mesh belonging to the function space you want to restrict.
Yes, in general, I would use locate_entities_boundary, and pass that information to locate_dofs_topological.

As you are asking a very generic question it is a bit hard to help you here.
if v1, x1, u1 all are based on the submesh, you should get no error with the code above (depending on how you defined function(x1, u1)).
This is because a submesh is just a mesh, and any operation on operators that only relate to it should work.

Please bear in mind that this post has diverged far from its original topic.
For further questions regarding the problem you describe above, please make a new post with a minimal reproducible example.