Add spring energy/condition

variational formulation
1

Hello,

I’m working under FEniCS for 1 year for my PhD thesis and i wanted to first of all thanks you for this amazing work.

I’m coding an hyper-elastic FE with a damage formulation and i want to add an energy (that does not require integration) to my minimization.

Lets define S1 the model domain, S2 the domain of a left suture and S3 for a right one.

When the model is stretched, some sutures need to be took into account in order to hold the sample in the perpendicular direction : Wwire = 0.5*kwire(u(x)dS2-u(x)dS3+l0)**2 with l0 the initial distance between suture.

I guess the code will be more clear that me :

# Function used to efficiently project CG1 on DG0 space (for H)
def local_project(v, V, u=None):
    dv = TrialFunction(V)
    v_ = TestFunction(V)
    a_proj = inner(dv, v_)*dx
    b_proj = inner(v, v_)*dx
    solver = LocalSolver(a_proj, b_proj)
    solver.factorize()
    if u is None:
        u = Function(V)
        solver.solve_local_rhs(u)
        return u
    else:
        solver.solve_local_rhs(u)
        return

### Load boundaries and subdomains defined under gmsh
subdomains = MeshFunction("size_t", mesh, path + "/mesh_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, path + "/mesh_facet_region.xml")
File(path + '/subdomains.pvd') << subdomains
File(path + '/boundaries.pvd') << boundaries

### Define domains
dx = Measure("dx", domain=mesh, subdomain_data=subdomains)
ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
dS = Measure("dS", domain=mesh, subdomain_data=boundaries)

sutleft.append(ds(8000))
sutright.append(ds(8100))
dstop = ds(8103)
dsright = ds(8104)
dsbottom = ds(8105)
dsleft = ds(8106)

### Define function space
Vu = VectorFunctionSpace(mesh, "CG", 1)
Vtens = TensorFunctionSpace(mesh, "DG", 0)
Vfunc = FunctionSpace(mesh, "DG", 0)

### Define Dirichlet boundary
bcbottom = Constant((0., 0.))
bctop = Constant((0., 0.))

bcs = [
        DirichletBC(Vu, bcbottom, boundaries, 8105),
        DirichletBC(Vu, bctop, boundaries, 8103)]

### Define functions
du = TrialFunction(Vu)            # Incremental displacement
v  = TestFunction(Vu)             # Test function

u = Function(Vu, name="Displacement")                
E = Function(Vtens, name="E")
S = Function(Vtens, name="PK1 Stress")

### Initial parameters
zeros = Constant((0.,  0.))
b = zeros  # Body force per unit volume
ttop = zeros  # Traction force on the boundary
tbottom = zeros  # Traction force on the boundary
tleft = zeros
tright = zeros

### Suture wire energy
kwire = Constant(400)
l0sut = Constant(3)
Wwire = 0.5*kwire*(u.sub(0)*sutleft-u.sub(0)*sutright + l0sut)

### Stored strain energy density
psi = PSI(u=u, ...) # computed with another script

### External volumetric force
Fvol = dot(b, u)*dx

### External surface force
Fsurf = dot(ttop, u)*dstop+dot(tbottom, u)*dsbottom+dot(tleft, u)*dsleft+dot(tright, u)*dsright

### Total potential energy
Pi = psi*dx - Fvol - Fsurf + Wwire

### Solve variational problem at intial conditions
u.interpolate(Constant((0., 0.)))
F = derivative(Pi, u, v)
J = derivative(F, u, du)

problem_u = NonlinearVariationalProblem(F, u, bcs=bcs, J=J)
solver_u = NonlinearVariationalSolver(problem_u)
solver_u.parameters['newton_solver']['linear_solver'] = 'gmres'
solver_u.parameters['newton_solver']['preconditioner'] = 'hypre_euclid'
solver_u.parameters['newton_solver']['relative_tolerance'] = 1e-4
solver_u.parameters['newton_solver']['maximum_iterations'] = 100
solver_u.parameters['newton_solver']['error_on_nonconvergence'] = False
solver_u.parameters['newton_solver']['krylov_solver']['relative_tolerance'] = 1e-4
solver_u.parameters['newton_solver']['krylov_solver']['nonzero_initial_guess'] = True
solver_u.parameters['symmetric'] = True

step = np.arange(0,1.1,0.1)
res_f = []
for i in step :
    bctop.assign(Constant((0., 30.*i)))
    solver_u.solve()

    ### Write solutions
    file_results.write(u, i)

    local_project(strain(u), Vtens, E)
    file_results.write(E, i)

    local_project(stress(u, SEF, x, Name), Vtens, S)
    file_results.write(S, i)

    local_project(psi, Vfunc, SSE)
    file_results.write(SSE, i)

I removed the damage part and i let the definition of the external energy to use it later.
As you can imagine i have some trouble to define Wwire and add it to the functional energy Pi.

i have read an old post a bit similar : https://fenicsproject.org/qa/6580/spring-functionals-assembling-non-integral-functional/ but without solution.

Thank you for the help