I guess the first question is, what is the reason for the four subdomains? If the objective is simply to have different properties (pde coefficients) in each subdomain, there is another way to do that. You can generate a single mesh with poimesh, selecting nx and ny so that you have element boundaries at the inner subdomain boundaries. Then you can use a function to define the coefficients in each element based on the x, y location of its centroid.
Bill
