# My Plane Strain Problem is not yielding any result

My problem of Plane Strain condition is not yielding any results even though it solves. I am not understanding if there is a problem with Plane Strain Formulation of some thing else.

``````import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import sys

from mpi4py import MPI
from petsc4py import PETSc
import pyvista

from dolfinx.cpp.mesh import to_type, cell_entity_type
#from dolfinx.fem import Constant, Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc
#from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector, create_vector, create_matrix, set_bc
#from dolfinx.fem.petsc import LinearProblem
from dolfinx import mesh, fem, io
import dolfinx.fem.petsc
from dolfinx.geometry import bb_tree, compute_collisions_points, compute_colliding_cells
from dolfinx.io import VTXWriter, distribute_entity_data, gmshio
from dolfinx.mesh import create_mesh, meshtags_from_entities
from dolfinx import default_scalar_type

from ufl import (FacetNormal, FiniteElement, as_matrix, Identity, Measure, TestFunction, TrialFunction, VectorElement, as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym)
import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem, MultiPointConstraint
``````
``````gmsh.initialize()

gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
h = 1

lc = 1e-3

# Added Points in the geometry

# Added Lines in the geometry

gmsh.model.geo.synchronize()

gmsh.option.setNumber("Mesh.CharacteristicLengthMin", h)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", h)

gmsh.model.mesh.generate(gdim)

domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
#facet_markers.name = "Facet markers"

gmsh.write("Winding Pack-1.msh")
gmsh.fltk.run()
gmsh.finalize()

``````
``````import ufl

def strain(u, repr ="vectorial"):
if repr =="vectorial":
return as_vector([eps_t[0,0], eps_t[1,1], 2*eps_t[0,1]])
elif repr =="tensorial":
return eps_t

E = fem.Constant(domain, 20.0)
nu = fem.Constant(domain,0.3)

lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

C = as_matrix([[lmbda+2*mu, lmbda, 0.0],[lmbda, lmbda+2*mu,  0.0],[0.0, 0.0, mu]])

def stress(u, repr ="vectorial"):
sigv = dot(C, strain(u))
if repr =="vectorial":
return sigv
elif repr =="tensorial":
return as_matrix([[sigv[0], sigv[2]], [sigv[2], sigv[1]]])

# Define Function Space
degree = 2
V = fem.functionspace(domain, ("P",degree, (gdim,)))

#Define Variational Problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = fem.Function(V, name = "Displacement")
a_form = inner(stress(du),strain(u_))*dx

# Applying a uniform pressure
#rho = 7850
#g = 9.81
T = fem.Constant(domain, 1000.0)
#fg = fem.Constant(domain, np.array([0, -rho*g]))

#Self-weight on the surface
n = FacetNormal(domain)
#L_form = dot(-fp*n,u_) * ds(5) + dot(fg*n,u_) * ds(5)

L_form = dot(T*n,u_) * ds(3)

#Boundary Conditions (simply supported on top and bottom edges)

V_ux, _ = V.sub(0).collapse()
V_uy, _ = V.sub(1).collapse()
bottom_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(1))
#top_dofs = fem.locate_dofs_topological((V.sub(1), V_uy), gdim-1, facet_markers.find(3))
left_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(4))
right_dofs = fem.locate_dofs_topological((V.sub(0), V_ux), gdim-1, facet_markers.find(2))

ux0 = fem.Function(V_ux)
uy0 = fem.Function(V_uy)
bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1)),fem.dirichletbc(ux0, left_dofs, V.sub(0)),fem.dirichletbc(ux0, right_dofs, V.sub(0))]
#bcs = [fem.dirichletbc(uy0, bottom_dofs, V.sub(1))]
#------------------------------------------------------------------------------------------------------------------------------------------

#Cyclic Symmetry Condition
#def PeriodicBoundary(x):
#   return np.isclose(x[0], 1, atol=atol)

#facets = locate_entities_boundary(mesh, mesh.topology.dim-1, PeriodicBoundary)
#arg_sort = np.argsort(facets)
#mt

#-------------------------------------------------------------------------------------------------------------------------------------------

#def epsilon(v):

#def sigma(v):
#   return lmbda*tr(epsilon(v))*Identity(dim)+ 2*mu*epsilon(v)

WPsolve = fem.petsc.LinearProblem(a_form, L_form, u=u, bcs=bcs)
WPsolve.solve()

V0 = fem.functionspace(domain, ("DG", 0, (3,)))
sig_exp = fem.Expression(stress(u), V0.element.interpolation_points())
sig = fem.Function(V0, name="Stress")
sig.interpolate(sig_exp)

vtk = io.VTKFile(domain.comm, "WPsolve_elasticity9.pvd", "w")
vtk.write_function(u, 0)
vtk.write_function(sig, 0)
vtk.close()
``````