# Result visualization for uniaxial loading of a cube with two different maerials

**URL:** <https://fenicsproject.discourse.group/t/result-visualization-for-uniaxial-loading-of-a-cube-with-two-different-maerials/5299>\
**Category:** Uncategorized\
**Created:** [March 13, 2021, 5:39pm UTC](https://fenicsproject.discourse.group/t/result-visualization-for-uniaxial-loading-of-a-cube-with-two-different-maerials/5299 "2021-03-13T17:39:00Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![swatisharmapanda](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/swatisharmapanda/32/109_2.png) [@swatisharmapanda](https://fenicsproject.discourse.group/u/swatisharmapanda)\
**Post date:** [March 13, 2021, 5:39pm UTC](https://fenicsproject.discourse.group/t/result-visualization-for-uniaxial-loading-of-a-cube-with-two-different-maerials/5299/1 "2021-03-13T17:39:00Z")

</div>

HI,  
I am trying to do uniaxial testing of a cube with two different nearly-incompressible NeoHookean materials for two sub-domains. My code is running fine. However, I do not know how to visualize the result. After the solver.solve() command, I calculate PK stress, but I can not figure out how to display both stresses in the mesh (As the mesh has two different materials, it should have two different stressed regions). Here is my code.

```
from dolfin import *
import matplotlib.pyplot as plt
import os
import dolfin as dlf
import numpy as np
import math 
from mshr import *

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"

# Kinematics
def pk1Stress(u,c1,nu):
    G = 2*c1
    K = 2*G*(1+nu)/(3*(1-2*nu)) # Bulk modulus
    I = Identity(V.mesh().geometry().dim()) # Identity tensor
    F = I + grad(u) # Deformation gradient ( I+ strain = L/L + dL/L = l/L = F)
    C = F.T*F # Right Cauchy-Green tensor
    Ic, J = tr(C), det(F) # Invariants of deformation tensors

    pk1 = 2*c1*F + (K*ln(J) - 2*c1) *inv(F).T # PK1 stress
    return pk1

def geometry_cube():
    mesh = UnitCubeMesh(4, 4, 4) 
    boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
    x0 = AutoSubDomain(lambda x: near(x[0], 0))
    x1 = AutoSubDomain(lambda x: near(x[0], 1))
    y0 = AutoSubDomain(lambda x: near(x[1], 0))
    z0 = AutoSubDomain(lambda x: near(x[2], 0))
    x0.mark(boundary_parts, 1)
    y0.mark(boundary_parts, 2)
    z0.mark(boundary_parts, 3)
    x1.mark(boundary_parts, 4)
    
    # Define a MeshFunction over two subdomains
    subdomains = MeshFunction('size_t', mesh, mesh.topology().dim())
    class Omega0(SubDomain):
        def inside(self, x, on_boundary):
            return True if x[1] <= 0.5 else False
    class Omega1(SubDomain):
        def inside(self, x, on_boundary):
            return True if x[1] >= 0.5 else False
        
    # Mark subdomains with numbers 0 and 1
    subdomain0 = Omega0()
    subdomain0.mark(subdomains, 0)
    subdomain1 = Omega1()
    subdomain1.mark(subdomains, 1)
    
    return boundary_parts,subdomains

# Create mesh and define function space ============================================
facet_function,subdomains = geometry_cube()
mesh = facet_function.mesh()
gdim = mesh.geometry().dim()
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
print('Number of nodes: ',mesh.num_vertices())
print('Number of cells: ',mesh.num_cells())

# Limit quadrature degree
dx = dx(degree=4)

# Create function space
V = VectorFunctionSpace(mesh, "CG", 1) 
du, _u = TrialFunction(V), TestFunction(V)
u = Function(V, name='u') 

# Create tensor function spaces for stress and stretch output 
W_DFnStress = TensorFunctionSpace(mesh, "CG", degree=1)
defGrad = Function(W_DFnStress, name='F')
PK1_stress = Function(W_DFnStress, name='PK1')

# Define Dirichlet boundary (symmetric Bc- no out of plane motion)
bc0 = DirichletBC(V.sub(0), Constant(0.), facet_function, 1)# ux = 0 at x =0
bc1 = DirichletBC(V.sub(1), Constant(0.), facet_function, 2)# uy = 0 at y =0
bc2 = DirichletBC(V.sub(2), Constant(0.), facet_function, 3)# uz = 0 at z =0
x1disp = Constant(0.0)
bc3 = DirichletBC(V.sub(0), x1disp, facet_function, 4)# prescribed disp 
bcs = [bc0,bc1,bc2,bc3]

# material parameters for different subdomains, NeoHookean material
c0 = 1000  
c1 = 2000
nu =0.49

# # Total potential energy
pkstrs0 = pk1Stress(u,c0,nu) #stress in sub domain 0
pkstrs1 = pk1Stress(u,c1,nu) #stress in sub domain 1

F1_0 = inner(pkstrs0, grad(_u))*dx(0) 
F1_1 = inner(pkstrs1, grad(_u))*dx(1) 

FF = F1_0+F1_1
DJ = derivative(FF,u,du)  

# Create nonlinear variational problem and solve
problem = NonlinearVariationalProblem(FF, u, bcs=bcs, J=DJ)
solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
solver.parameters['newton_solver']['linear_solver'] = 'mumps'

# Time stepping parameters
Tend = 1
Nsteps = 10
time = np.linspace(0, Tend, Nsteps+1)
dt = time[1]-time[0]

# Save solution in xdmf format for visualization using paraview
file_results = XDMFFile("./Results/TestCubeStretching/Uniaxial.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

for i in range(len(time)):
    tt = time[i]
    print(i, 'time: ', tt)

    # Increase traction
    x1disp.assign(0.5*tt)

    # solve and save disp
    solver.solve()

    # Save DefGrad to file in xdmf format
    I = Identity(V.mesh().geometry().dim()) # Identity tensor
    DF = I + grad(u)
    defGrad.assign(project(DF, W_DFnStress))

    # Save Stress to file in VTK format (I have an issue here. for stress calculation what expression should I use)
    PK1t= pk1Stress(u,c0,nu) # PK1t= pk1Stress(u,c1,nu) which expression to use ????
    PK1_stress.assign(project(PK1t, W_DFnStress))

    # save xdmf file
    file_results.write(u, tt)
    file_results.write(defGrad, tt)
    file_results.write(PK1_stress, tt)

print('done')
```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [March 14, 2021, 4:50pm UTC](https://fenicsproject.discourse.group/t/result-visualization-for-uniaxial-loading-of-a-cube-with-two-different-maerials/5299/2 "2021-03-14T16:50:22Z")

</div>

You are projecting your stress into a CG 1 space, when the stressed are contained in a DG0 space. You should project them into an appropriate function space and use `XDMFFile.write_checkpoint` to write it to a suitable file format for visualization in paraview. See for instance: [When using CG degree 1, the result corresponds to degree 2 - #10 by dokken](https://fenicsproject.discourse.group/t/when-using-cg-degree-1-the-result-corresponds-to-degree-2/2751/10)