Stress Distribution over 2 material mesh

General
1

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
      1 #grad_u= project(eps_t(v)[1],W)
      2 
      3 #plot(grad_u[1])
----> 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 ?
    Please find the MWE.
def str(v):
    E1= as_vector([0,v[1].dx(0),v[1].dx(0)])
    return E1

grad_u= project(str(w),W)
plot(grad_u[1])

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?