Error implementing a constraint on a 2D crack propagation problem

Linear Algebra
1

Hello,

I am trying to solve the 2-D elasticity equations coupled with the phase field equations to look into crack propagation. I have 2 variational equations : one for the displacement and the other for the phase field. The idea is that cracking is allowed in tension but not in compression. I am running a loop on the number of vertices to do the following :

  1. Solving the displacement problem to get the new displacement.
  2. Updating the maximum positive strain energy
  3. Solving the phase field equation with the updated values
  4. Implement a constraint as shown below:

In step no 4. I used the following:

Psi0P_TS = LocalProject(Psi0P(unew, Q), Q)
Psi0N_TS = LocalProject(Psi0N(unew, Q), Q)

Psi0P_TS_Array = Psi0P_TS.vector().vec().getArray()
Psi0N_TS_Array = Psi0N_TS.vector().vec().getArray()

m_Array = m.vector().vec().getArray()

for j in range( mesh.num_vertices() ):

	if ( Psi0P_TS_Array[j] < Psi0N_TS_Array[j] ):
	
		m_Array[j] = 0.

m.vector().set_local(m_Array)

assign(mold,m)

LocalProject is just a user defined function that gives me the projected values on the functionspace.
The idea is that it should check at every point whether the positive strain energy is greater than the negative strain energy - if not, the phase field parameter(m) is reduced to zero.

On running this, however, it gives me an error. Either the solution does not converge or even if it is made to converge using low tolerance values - the result shows a static crack. On the other hand, if I remove this constraint - it shows the crack propagating (even though it is wrong) .

Please help me out. I can’t see where I am going wrong. Thanks in advance!