Two materials and Mechanism

General
4

Its not very clear what you are asking.
You are not using dx1 and dx2 in the correct way.

Note that your example is far from minimal (as you could have stopped the whole example after generating vm, and visualize it in Paraview to see if it is correct.

Here is a complete example, that also includes the mesh conversion.

from fenics import *
import meshio
from pathlib import Path
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

mesh_file = Path("mecanismo_mesh.msh")
xdmf_file = mesh_file.with_suffix(".xdmf")
facet_file = xdmf_file.with_stem("mecanismo_facet")

msh = meshio.read(mesh_file)
mesh3D = create_mesh(msh, "tetra")
mesh2D = create_mesh(msh, "triangle")
meshio.write(xdmf_file, mesh3D)
meshio.write(facet_file, mesh2D)


# Load
F = 2e2

# Define 3D geometry
mesh = Mesh()
with XDMFFile(str(xdmf_file)) as infile:
    infile.read(mesh)  
mvc2 = MeshValueCollection("size_t", mesh, 2) 
with XDMFFile(str(facet_file)) as infile:
    infile.read(mvc2, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
mvc3 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile(str(xdmf_file)) as infile:
    infile.read(mvc3, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc3)


# Define Material
Em = [210e3, 100]
#Em = [100, 210e3] # inverting

Vm = FunctionSpace(mesh, "DG", 0)
vm = Function(Vm)

for cell in cells(mesh):
    i = cell.index()
    if cf.array()[i] == 22:
        vm.vector().vec().setValueLocal(i, Em[0])
    elif cf.array()[i]==23:
        vm.vector().vec().setValueLocal(i, Em[1])

# Define function space and base functions
V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)

# Boundary Condition
bc = DirichletBC(V, Constant((0.0, 0.0, 0.0)), mf, 20)

BCs = [bc]

# Load
ds = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=21)
t = Expression(('q_x', 'q_y', 'q_z'), degree = 1, q_x = 0.0, q_y = 0.0, q_z = -F)

## Constitutive relation

def eps(v):
    return sym(grad(v))

#E = Constant(210e3)
E = vm
nu = Constant(0.3)
model = "plane_stress"
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
    lmbda = 2*mu*lmbda/(lmbda+2*mu)

def sigma(v):
    dim = v.geometric_dimension()
    return lmbda*tr(eps(v))*Identity(dim) + 2.0*mu*eps(v)

## Variational formulation

du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
L = dot(t, u_)*ds

A = assemble(a)
B = assemble(L)

[bc.apply(A, B) for bc in BCs] 

u = Function(V, name="Displacement")
solve(A ,u.vector() ,B)

#file_results
file_results = XDMFFile("u_two_material.xdmf")
#file_results = XDMFFile("u_two_material_inverting.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
file_results.close()

Choosing the inverse values gives a significantly different solution.

For any further inquiries, please make an effort in isolating what you have problems with (in this case, the definition of E through vm).