# Submesh / FunctionSpace for subdomain

**URL:** <https://fenicsproject.discourse.group/t/submesh-functionspace-for-subdomain/11287>\
**Category:** mesh\
**Created:** [May 25, 2023, 9:12pm UTC](https://fenicsproject.discourse.group/t/submesh-functionspace-for-subdomain/11287 "2023-05-25T21:12:02Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\
**Post date:** [May 25, 2023, 9:12pm UTC](https://fenicsproject.discourse.group/t/submesh-functionspace-for-subdomain/11287/1 "2023-05-25T21:12:02Z")

</div>

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](https://comet-fenics.readthedocs.io/en/latest/demo/periodic_homog_elas/periodic_homog_elas.html)

```auto
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.

```auto
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

```auto
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!

---

<div class="post-metadata">

**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\
**Post date:** [May 25, 2023, 11:23pm UTC](https://fenicsproject.discourse.group/t/submesh-functionspace-for-subdomain/11287/2 "2023-05-25T23:23:38Z")

</div>

**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)?

```auto
# 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

> [@Usage of FunctionAssigner](https://fenicsproject.discourse.group/t/usage-of-functionassigner/4868/1):
>
> ```auto
> T = FunctionSpace(mesh, 'CG', 2)
> 
> v = Function(V)
> 
> S_1 = Function(S)
> S_11, S_12 = S_1.split()
> 
> T_S_1 = interpolate(S_11, T)
> T_S_2 = interpolate(S_12, T)
> #next line results in an error
> assigner = FunctionAssigner(V, [T, T])
> assigner.assign(v, [T_S_1, T_S_2])
> 
> ```

---

<div class="post-metadata">

**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\
**Post date:** [May 26, 2023, 2:17pm UTC](https://fenicsproject.discourse.group/t/submesh-functionspace-for-subdomain/11287/3 "2023-05-26T14:17:19Z")

</div>

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

> [@bagla0](#):
>
> Periodic homogenization of linear elastic materials — Numerical tours of continuum mechanics using FEniCS master documentation

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

> [@bagla0](#):
>
> ```auto
> 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
> 
> ```

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.

> [@bagla0](#):
>
> ```auto
> # 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)
> 
> ```

The plots obtained are

> [@bagla0](#):
>
> ```auto
> plt.figure()
> p1 = plot(sigm)
> plt.colorbar(p1)
> plt.show()
> 
> ```

![mat](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/d/d9b9273d4b0a845e4ffcb8dfdc8c763dea50a2b8.jpeg)

> [@bagla0](#):
>
> ```auto
> plt.figure()
> p2 = plot(sigf)
> plt.colorbar(p2)
> plt.show()
> 
> ```

![fib](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/e/e98ce692f55a02558255005711003e0567e015eb.png)

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