Submesh / FunctionSpace for subdomain

mesh
1

Hello FEniCS family,
Problem in finding Functionspace specific to subdomain region of mesh to project sigma.

For periodic homogenization tutorial, Periodic homogenization of linear elastic materials — Numerical tours of continuum mechanics using FEniCS master documentation

dx = Measure('dx')(domain=mesh, subdomain_data=subdomains)

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,[])

def sigma_vec(v, i,Eps):
        
     E,nu=material_parameters[i]     
     if i==0:  
        lmbda = E*nu/((1+nu)*(1-2*nu))
        mu = E/(2*(1+nu))
        C1=lmbda+2*mu
        C=as_tensor([(C1,lmbda,lmbda,0,0,0),(lmbda,C1,lmbda,0,0,0),(lmbda,lmbda,C1,0,0,0),(0,0,0,mu,0,0),(0,0,0,0,mu,0),(0,0,0,0,0,mu)])
        s1= dot(C,eps(v)+Eps)
    else:
        lmbda = E*nu/((1+nu)*(1-2*nu))
        mu = E/(2*(1+nu))
        C1=lmbda+2*mu
        C=as_tensor([(C1,lmbda,lmbda,0,0,0),(lmbda,C1,lmbda,0,0,0),(lmbda,lmbda,C1,0,0,0),(0,0,0,mu,0,0),(0,0,0,0,mu,0),(0,0,0,0,0,mu)])
        s1= dot(C,eps(v)+Eps)
    return s1

I want to find the stress distribution through project command.

sigma_w = project(sigma_vec(w,0,Eps), W)

where 0 defines subdomain 0. Its calculated sigma using only one material stiffness matrix.

I tried with

sigma_w = project(sum([sigma_vec(w,i,Eps) for i in range(nphases)]), W)

but, this gave summation of stress at each element taken both material properties seperately.
How can I get the sigma projection taking respective material stiffness matrix based on its subdomain of 0 or 1 ?

I realize, this coul be done using Submesh feature to define 2 different Function Space, but, didn’t know exactly what to do.
Any help is greatly appreciated!

2

UPDATE

I am able to generate the submesh and get stress plot for each submesh domain seperately.
How can I merge these plots (using hold on feature similar to matlab)?

# stress Plot
Vf = FunctionSpace(submesh_fib, 'DG', 0)
Vm = FunctionSpace(submesh_mat, 'DG', 0)
sigf=project(sigma_vec(v,1,Eps)[0],Vf)
sigm=project(sigma_vec(v,0,Eps)[0],Vm)

plt.figure()
p1 = plot(sigm)
plt.colorbar(p1)
plt.show()


plt.figure()
p2 = plot(sigf)
plt.colorbar(p2)
plt.show()
![fib|569x451](upload://xk57ZBprNVpqXI9bTeojdL3vrPd.png)

I tried FunctionAssigner, but getting error in interpolate step according to

3

I am stuck at this point. I have generated two submesh region based on periodic homogenization tutorial 2D RVE.

The stress function below takes fluctuation function, v (function output of linear variational solver).

Two function space were crated based on two submesh and stress were plotted in each submesh function space.
I want to get a single plot for two function space.

The plots obtained are

mat

fib

Kindly help me to get the common plot. I apologize if I am unable to explain my query.