Error in results for spatially varying material property

[abhiawasthi1993]

Hello, I am attempting to solve a 2D plane strain boundary value problem using a linear elastic material with Young's modulus that varies spatially, expressed as E = E(x,y) or E(nodes). I have used 'QUAD4' as element type. The basis functions utilized for the solution are assumed to support the Young's modulus. Below is the code I have employed to convert the values from nodes to the elements:

class Elasticity(FEM):
    def custom_init(self):
        """ Override base class method
        """
        self.flex_inds = np.arange(len(self.points))
        
    def get_tensor_map(self):
        def stress(u_grad, theta):
            nu = 0.3
            E = theta[0]
            epsilon = 0.5*(u_grad + u_grad.T)
            eps11 = epsilon[0,0]
            eps22 = epsilon[1,1]
            eps12 = epsilon[0,1]
            sig11 = E/((1 + nu)*(1 - 2*nu))*((1-nu)*eps11 + nu*eps22)
            sig22 = E/((1 + nu)*(1 - 2*nu))*(nu*eps11 + (1-nu)*eps22) 
            sig12 = E/((1 + nu)*2)*eps12
            sigma = np.array([[sig11, sig12],[sig12, sig22]])
            return sigma
        return stress
    
    def set_params(self, params):
        NCA = self.cells
        NCA = np.array(NCA)
        
        def compute_params(e):
            a = NCA[e,:]
            Ee = params[a]
            E = np.matmul(self.shape_vals, Ee)
            return np.transpose(E)
        
        result = jax.vmap(compute_params)(np.arange(self.num_cells))
        result = np.transpose(result, axes=(0,2,1))
        
        self.full_params = params
        self.internal_vars['laplace'] = [result]

When I used a uniform material parameter distribution (param = 0.01*np.ones((len(problem.flex_inds), 1))), the results matched those from a benchmark study. However, when I introduced a spatially varying (heterogeneous) distribution for the material parameter, the accuracy of the results decreased.

I attempted to find solutions in existing examples, but most of them assumed that the varying field had a constant value within each element, which does not apply to my case.

I would appreciate any guidance or suggestions to identify where I might have made errors in my implementation for handling the spatially varying material parameter and to improve the accuracy of the results.