How to calculate two max values of vector in two different subdomains?

General
1

Hi, I have unit square as my domain and divided it into two subdomains with different material properties using MeshFunction(), I want to calculate the maximum value of \sigma_{yy} (sigma(u)[1, 1] in my case) differently for different subdomains having 2 different material properties. I am able to calculate the maximum value over whole domain using :-

stress = sigma(u)[1, 1]
Vs = FunctionSpace(mesh, 'DG', 0)
stress = local_project(stress, Vs)
print(stress.vector().max())

But not able to do it for 2 subdomains separately.
Here is my MWE :-

from fenics import *
from mshr import *
import numpy as np

domain = Rectangle(Point(0, 0),Point(1, 1))
mesh = generate_mesh(domain, 10)

V = VectorFunctionSpace(mesh, 'CG', 2)

tol = 1E-14
class Omega_0(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] <= 0.5 + tol

class Omega_1(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] >= 0.5 - tol

materials = MeshFunction('size_t', mesh, 2)

subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(materials, 0)
subdomain_1.mark(materials, 1)

materials.set_all(0)
subdomain_1.mark(materials, 1)

class al(UserExpression):
    def __init__(self, materials, al0, al1, **kwargs):
        super().__init__(**kwargs)
        self.materials = materials
        self.k_0 = al0
        self.k_1 = al1
    def eval_cell(self, values, x, cell):
        if self.materials[cell.index] == 0:
            values[0] = self.k_0
        else:
            values[0] = self.k_1
    def value_shape(self):
        return ()      
    
E1 = 21e3 
nu1 = 0.33
E2 = 210e3 
nu2 = 0.25

E = al(materials, E1, E2, degree = 0)
nu = al(materials, nu1, nu2, degree = 0)

lambda_ = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)

class bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0., 1e-8)
    
class top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 1., 1e-8)

boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)

bottom().mark(boundaries, 1)
top().mark(boundaries, 2)

bc_top = DirichletBC(V.sub(1), Constant(0.01), boundaries, 2)
bc_bottom = DirichletBC(V.sub(1), Constant(0.), boundaries, 1)

bc = [bc_bottom, bc_top]

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)

def sigma(u):
    return lambda_*tr(epsilon(u))*Identity(d) + 2*mu*epsilon(u)

u = TrialFunction(V)
d = u.geometric_dimension()
v = TestFunction(V)
f = Constant((0, 0))
T = Constant((0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

def local_project(v, V, u=None):
    dv = TrialFunction(V)
    v_ = TestFunction(V)
    a_proj = inner(dv, v_)*dx
    b_proj = inner(v, v_)*dx
    solver = LocalSolver(a_proj, b_proj)
    solver.factorize()
    if u is None:
        u = Function(V)
        solver.solve_local_rhs(u)
        return u
    else:
        solver.solve_local_rhs(u)
        return
    
u = Function(V)
solve(a == L, u, bc)

stress = sigma(u)[1, 1]
Vs = FunctionSpace(mesh, 'DG', 0)
stress = local_project(stress, Vs)
print(stress.vector().max())

I would appreciate any pointers in the direction of how I could achieve 2 different max values in 2 different subdomains.
Thanks.

2

You could do something along the lines of:

one_indices = materials.where_equal(1)
zero_indices = materials.where_equal(0)
print(stress.vector().get_local()[one_indices].max())
print(stress.vector().get_local()[zero_indices].max())

Note that to do this in parallel, you would have to communicate the max value of each process to find a global one.

3

This is exactly what I was looking for. Thanks.