My Plane Strain Problem is not yielding any result

dolfinx
1

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.graph import adjacencylist
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()

gmsh.model.add("Winding Pack-1")

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

lc = 1e-3

# Added Points in the geometry

gmsh.model.geo.addPoint(0,0,0,lc,1)
gmsh.model.geo.addPoint(1,0,0,lc,2)
gmsh.model.geo.addPoint(1.5,1,0,lc,3)
gmsh.model.geo.addPoint(-0.5,1,0,lc,4)
#gmsh.model.geo.addPoint(0.3,0.3,0,lc,5)
#gmsh.model.geo.addPoint(0.7,0.3,0,lc,6)
#gmsh.model.geo.addPoint(0.7,0.6,0,lc,7)
#gmsh.model.geo.addPoint(0.3,0.6,0,lc,8)

# Added Lines in the geometry

gmsh.model.geo.addLine(1,2,1)
gmsh.model.geo.addLine(2,3,2)
gmsh.model.geo.addLine(3,4,3)
gmsh.model.geo.addLine(4,1,4)
#gmsh.model.geo.addLine(5,6,5)
#gmsh.model.geo.addLine(6,7,6)
#gmsh.model.geo.addLine(7,8,7)
#gmsh.model.geo.addLine(8,5,8)

gmsh.model.geo.addCurveLoop([1,2,3,4],1)

gmsh.model.geo.addPlaneSurface([1],1)

gmsh.model.geo.synchronize()

gmsh.model.addPhysicalGroup(1, [1], name="bottom")
gmsh.model.addPhysicalGroup(1, [2], name="right")
gmsh.model.addPhysicalGroup(1, [3], name="top")
gmsh.model.addPhysicalGroup(1, [4], name="left")

gmsh.model.addPhysicalGroup(2, [1], name="my_surface")

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"):
    eps_t = sym(grad(u))
    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):
 #   return sym(grad(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()