# Stress Distribution over 2 material mesh

Hi FEniCS family,
My work is related to periodic homogenization example Periodic homogenization of linear elastic materials — Numerical tours of continuum mechanics using FEniCS master documentation.

I want to get the stress and strain distribution over RVE after getting my fluctuation functions, w through solve command. The following code is taken from that example in above link.

``````Ve = VectorElement("CG", mesh.ufl_cell(), 2)
Re = VectorElement("R", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10))
V = FunctionSpace(mesh, Ve)

v_,lamb_ = TestFunctions(W)
dv, dlamb = TrialFunctions(W)
w = Function(W)
dx = Measure('dx')(subdomain_data=subdomains)

Eps = Constant(((1, 0), (0, 0)))
F = sum([inner(sigma(dv, i, Eps), eps(v_))*dx(i) for i in range(nphases)])
a, L = lhs(F), rhs(F)
a += dot(lamb_,dv)*dx + dot(dlamb,v_)*dx
solve(a == L, w,[])
(v, lamb) = split(w)
``````
1. How can we go for finding stress distribution when we have two material phases?
(Like previous use sum during variational form define works, but how to use here)
I tried for single phase but, got the error in that too
``````stress_u= project(sigma(v,i,Eps)*dx(0),W)

Can only integrate scalar expressions. The integrand is a tensor expression with value shape (6,) and free indices with labels ().
---------------------------------------------------------------------------
UFLException                              Traceback (most recent call last)
Cell In[238], line 5
2
----> 5 stress_u= project(sigma(v,0,Eps)*dx(0),W)

File /usr/lib/python3/dist-packages/ufl/measure.py:408, in Measure.__rmul__(self, integrand)
406 # Allow only scalar integrands
407 if not is_true_ufl_scalar(integrand):
--> 408     error("Can only integrate scalar expressions. The integrand is a "
409           "tensor expression with value shape %s and free indices with labels %s." %
410           (integrand.ufl_shape, integrand.ufl_free_indices))
412 # If we have a tuple of domain ids build the integrals one by
413 # one and construct as a Form in one go.
414 subdomain_id = self.subdomain_id()

File /usr/lib/python3/dist-packages/ufl/log.py:158, in Logger.error(self, *message)
156 "Write error message and raise an exception."
157 self._log.error(*message)
--> 158 raise self._exception_type(self._format_raw(*message))

UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (6,) and free indices with labels ().
``````
1. For strain distribution, I used the project command which is working to give colored plot but incorrect. How can we get the correct strain distribution plot ?
``````def str(v):
E1= as_vector([0,v[1].dx(0),v[1].dx(0)])
return E1

Object cannot be plotted directly, projecting to piecewise linears.
<matplotlib.tri._tricontour.TriContourSet at 0x7f5dc70d0fa0>
``````
1. After split command ` (v,lam)=split(w)` , w is a coefficient whose nodal values can be computed using ` w.compute_nodal_values()` but, v and str(v) have becomes list tensor whose nodal values or values at guass point cannot be determined. Do we need to convert v to Coefficient?

2. I have seen in few post processing examples where Function space is defined with ‘P’ or ‘DG’

``````Q = FunctionSpace(mesh, "DG", 0)
u_DG = project(w, Q)

Shapes do not match: <Argument id=140040738440576> and <Coefficient id=140040744598336>.
---------------------------------------------------------------------------
UFLException                              Traceback (most recent call last)
Cell In[78], line 2
1 Q = FunctionSpace(mesh, "DG", 0)
----> 2 u_DG = project(w, Q)

File /usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/projection.py:129, in project(v, V, bcs, mesh, function, solver_type, preconditioner_type, form_compiler_parameters)
127 Pv = TrialFunction(V)
128 a = ufl.inner(w, Pv) * dx
--> 129 L = ufl.inner(w, v) * dx
131 # Assemble linear system
132 A, b = assemble_system(a, L, bcs=bcs,
133                        form_compiler_parameters=form_compiler_parameters)

File /usr/lib/python3/dist-packages/ufl/operators.py:158, in inner(a, b)
156 if a.ufl_shape == () and b.ufl_shape == ():
157     return a * Conj(b)
--> 158 return Inner(a, b)

File /usr/lib/python3/dist-packages/ufl/tensoralgebra.py:147, in Inner.__new__(cls, a, b)
145 ash, bsh = a.ufl_shape, b.ufl_shape
146 if ash != bsh:
--> 147     error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)))
149 # Simplification
150 if isinstance(a, Zero) or isinstance(b, Zero):

File /usr/lib/python3/dist-packages/ufl/log.py:158, in Logger.error(self, *message)
156 "Write error message and raise an exception."
157 self._log.error(*message)
--> 158 raise self._exception_type(self._format_raw(*message))

UFLException: Shapes do not match: <Argument id=140040738440576> and <Coefficient id=140040744598336>.

``````

Why we choose either instead of our original defined function space? I am not finding any resource to resolve this.

1. How do we get stress distribution over guass points and nodal points seperately ?
I realize project command gives only over guass points?