Let’s say my mesh has 4 cells,
cell 1 and 2 have material A and
cell 3 and 4 have material B.
So i make 2 physical groups in gmsh to use as sub-domains in fenicsx.
but for post processing i want to view results only for cell 1 and 3 . how do i do that??
I tried making another physical group in gmsh comprising of cells 1 and 3. but when i read that mesh in fenicsx, kernel crashes.
I think fenicsx doesnt allow cells belonging to multiple sub-domains. Is there a work around this issue??
Use dolfinx.mesh.create_submesh to create a mesh comprising of only cells marked with 1 and 3, and interpolate the solution onto that mesh using: Interpolation to and from sub-meshes for co-dimension 0 by jorgensd · Pull Request #3114 · FEniCS/dolfinx · GitHub
import gmsh
import sys
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
import matplotlib.pyplot as plt
import dolfinx
gmsh.initialize()
gmsh.option.setNumber('General.Verbosity', 0)
gmsh.model.add("square")
gmsh.model.occ.addRectangle(0,0,0,1,1,1)
gmsh.model.occ.addRectangle(1,0,0,1,1,2)
gmsh.model.occ.addRectangle(0,1,0,1,1,3)
gmsh.model.occ.addRectangle(1,1,0,1,1,4)
gmsh.model.occ.fragment([(2,1)],[(2,2),(2,3),(2,4)],removeObject=True,removeTool=True)
gmsh.model.occ.synchronize()
lines = gmsh.model.getEntities(1)
surf = gmsh.model.getEntities(2)
for i in lines:
gmsh.model.mesh.setTransfiniteCurve(i[1],2)
for i in surf:
gmsh.model.mesh.setTransfiniteSurface(i[1])
gmsh.model.mesh.setRecombine(i[0],i[1])
gmsh.model.addPhysicalGroup(2,[1,2],101)
gmsh.model.setPhysicalName(2,101,"matA")
gmsh.model.addPhysicalGroup(2,[3,4],102)
gmsh.model.setPhysicalName(2,102,"matB")
# gmsh.model.addPhysicalGroup(2,[1,3],103)
# gmsh.model.setPhysicalName(2,103,"post_process")
gmsh.model.mesh.generate(2)
# if '-nopopup' not in sys.argv:
# gmsh.fltk.run()
domain, markers, facets = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)
V = fem.functionspace(domain, ("Lagrange", 1))
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)
Q = fem.functionspace(domain, ("DG", 0))
kappa = fem.Function(Q)
mA_cells = markers.find(101)
kappa.x.array[mA_cells] = np.full_like(mA_cells, 1, dtype=default_scalar_type)
mB_cells = markers.find(102)
kappa.x.array[mB_cells] = np.full_like(mB_cells, 0.1, dtype=default_scalar_type)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type(-6))
a = ufl.inner(kappa * ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["u"] = uh.x.array.real
u_grid.set_active_scalars("u")
u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=True)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
u_plotter.show()
so in this MWE,
tag 101 has cells 1,2 with material A,
tag 102 has cells 3,4 with material B
Doing this
submesh, entity_map, vertex_map, geom_map = mesh.create_submesh(domain, 2, markers.find(101))
will give me results for cells 1 and 2 right??
how do I make a submesh for cells 1 and 3 only??
Create a new meshtag object consisting of cell 1 and 3. hardcoded that would look like:
indices = np.array([1,3], dtype=np.int32)
values = np.full_like(indices, 1)
combined_tag = dolfinx.mesh.meshtags(domain, domain.topology.dim, indices, values)
submesh, entity_map, vertex_map, geom_map = mesh.create_submesh(domain, domain.topology.dim,combined_tag.find(1))
Thanks this worked!
How can i create sub_mesh using entity tags from gmsh?
Because, if I refine the above mesh, I no longer know the element ids.
I tried getting element ids from gmsh.model.mesh.getElements but the result was not what I expected.
here is the full code,
import gmsh
import sys
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
import matplotlib.pyplot as plt
import dolfinx
gmsh.initialize()
gmsh.option.setNumber('General.Verbosity', 0)
gmsh.model.add("square")
gmsh.model.occ.addRectangle(0,0,0,1,1,1)
gmsh.model.occ.addRectangle(1,0,0,1,1,2)
gmsh.model.occ.addRectangle(0,1,0,1,1,3)
gmsh.model.occ.addRectangle(1,1,0,1,1,4)
gmsh.model.occ.fragment([(2,1)],[(2,2),(2,3),(2,4)],removeObject=True,removeTool=True)
gmsh.model.occ.synchronize()
lines = gmsh.model.getEntities(1)
surf = gmsh.model.getEntities(2)
for i in lines:
gmsh.model.mesh.setTransfiniteCurve(i[1],5)
for i in surf:
gmsh.model.mesh.setTransfiniteSurface(i[1])
gmsh.model.mesh.setRecombine(i[0],i[1])
gmsh.model.addPhysicalGroup(2,[1,2],101)
gmsh.model.setPhysicalName(2,101,"matA")
gmsh.model.addPhysicalGroup(2,[3,4],102)
gmsh.model.setPhysicalName(2,102,"matB")
# gmsh.model.addPhysicalGroup(2,[1,3],103)
# gmsh.model.setPhysicalName(2,103,"post_process")
gmsh.model.mesh.generate(2)
elements = list(gmsh.model.mesh.getElements(2,1)[1][0]) + list(gmsh.model.mesh.getElements(2,4)[1][0])
# if '-nopopup' not in sys.argv:
# gmsh.fltk.run()
domain, markers, facets = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)
V = fem.functionspace(domain, ("Lagrange", 1))
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
# Create facet to cell connectivity required to determine boundary facets
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)
Q = fem.functionspace(domain, ("DG", 0))
kappa = fem.Function(Q)
mA_cells = markers.find(101)
kappa.x.array[mA_cells] = np.full_like(mA_cells, 1, dtype=default_scalar_type)
mB_cells = markers.find(102)
kappa.x.array[mB_cells] = np.full_like(mB_cells, 0.1, dtype=default_scalar_type)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type(-6))
a = ufl.inner(kappa * ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["u"] = uh.x.array.real
u_grid.set_active_scalars("u")
u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=True)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
u_plotter.show()
indices = np.array(elements, dtype=np.int32)
values = np.full_like(indices, 1)
combined_tag = dolfinx.mesh.meshtags(domain, domain.topology.dim, indices, values)
submesh, entity_map, vertex_map, geom_map = mesh.create_submesh(domain, domain.topology.dim,combined_tag.find(1))
Vs = fem.functionspace(submesh, ("Lagrange", 1))
uDs = fem.Function(Vs)
uDs.interpolate(uh,None,entity_map)
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(Vs)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["u"] = uDs.x.array.real
u_grid.set_active_scalars("u")
u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=True)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
u_plotter.show()
This is what i get in 4 element case (as desired)
This is what I get in refined mesh
Instead of doing this, i.e. creating two groups, I would advice you to create one group for each domain, i.e. 101->1, 102->2, 103->3, 104->4 etc and then combine them in DOLFINx rather than in GMSH (as it is expensive to sort out duplicate elements from gmsh meshes).
Say that you have tagged your domains as sketched above, you can easily set values for each of them, and then create a combined meshtags object for some of them:
combined_indices = []
for tag in [1,3]:
combined_indices.append(mt.find(tag))
combined_indices = np.sort(np.hstack(combined_indices))
submesh, entity_map, vertex_map, geom_map = mesh.create_submesh(domain, domain.topology.dim, combined_indices)
This works!
But having separate physical group for each domain leads to problems in large meshes.
Consider a 3D linear anisotropic elasticity,
for n domains I have n 6by6 stiffness matrices
So when I do this,
for i in range(n):
a += ufl.inner(ortho_sigma(u,C[i]), ortho_epsilon(v))*dx(i)
I start getting errors like
libffcx_forms_4483306f7a1bb6e036a016a5df08e04a32ac95a2.c:8399751: note: adding ‘-flarge-source-files’ will allow for more column-tracking support, at the expense of compilation time and memory
Which is why I was combining multiple domains with same stiffness matrix in a single physical group.
I have never seen this error before.
Please provide an MWE, as there might be way easier ways of doing what you would like to do.
I am not able to reproduce this error anymore.
But is this the best way to assign different material properties(stiffness matrix) to multiple domains
I would need to see the definition of C[i]. Please make an illustrative code.
Here’s the full code: Here C is constant for all subdomains for simplicity (Ideally C would vary as per subdomains)
import gmsh
import sys
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
import matplotlib.pyplot as plt
import dolfinx
import time
def make_layers(n):
length = 1000;width = 1000; thickness= 0.5
gmsh.initialize()
gmsh.option.setNumber("General.Verbosity",0)
gmsh.model.add("plate")
for i in range(n):
gmsh.model.occ.addBox(0,0,i*thickness,length,width,thickness)
gmsh.model.occ.synchronize()
vols = gmsh.model.getEntities(3)
for i in range(n-1):
gmsh.model.occ.fragment([vols[i]],[vols[i+1]],removeObject=True,removeTool=True)
gmsh.model.occ.synchronize()
lines = gmsh.model.getEntities(1)
lengths = []
for i in lines:
lengths.append(gmsh.model.occ.getMass(1,i[1]))
surf = gmsh.model.getEntities(2)
vols = gmsh.model.getEntities(3)
for i in range(len(lines)):
if lengths[i] == length:
gmsh.model.mesh.setTransfiniteCurve(lines[i][1],2)
else:
gmsh.model.mesh.setTransfiniteCurve(lines[i][1],0)
for i in surf:
gmsh.model.mesh.setTransfiniteSurface(i[1])
gmsh.model.mesh.setRecombine(i[0],i[1])
for i in vols:
gmsh.model.mesh.setTransfiniteVolume(i[1])
f_surfaces = gmsh.model.getEntitiesInBoundingBox(-length/2,-width/2,-thickness/2,length/2,(n+1)*width,(n+1)*thickness,2)
l_surfaces = gmsh.model.getEntitiesInBoundingBox(length/2,-width/2,-thickness/2,2*length,(n+1)*width,(n+1)*thickness,2)
fixed_surfaces = []
load_surfaces = []
vol_count = 1000
for i in range(n):
fixed_surfaces.append(f_surfaces[i][1])
load_surfaces.append(l_surfaces[i][1])
gmsh.model.addPhysicalGroup(3,[i+1],vol_count)
gmsh.model.setPhysicalName(3,vol_count,"volume"+ str(i+1))
vol_count += 1
gmsh.model.addPhysicalGroup(2,fixed_surfaces,101)
gmsh.model.setPhysicalName(2,101,"fixed_surface")
#print(fixed_surfaces)
gmsh.model.addPhysicalGroup(2,load_surfaces,102)
gmsh.model.setPhysicalName(2,102,"force_surface")
#print(load_surfaces)
gmsh.model.mesh.generate(3)
gmsh.model.mesh.refine()
# if '-nopopup' not in sys.argv:
# gmsh.fltk.run()
domain, markers, facets = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)
gmsh.finalize()
return domain,markers,facets
def fe_case(n):
domain,markers,facets = make_layers(n)
E = 1.1e9;nu = 0.3
mu = E/2/(1+nu);lambda_ = E*nu/(1+nu)/(1-2*nu)
C = ufl.as_matrix([[float((E/((1+nu)*(1-2*nu)))*(1-nu)), float((E/((1+nu)*(1-2*nu)))*nu), float((E/((1+nu)*(1-2*nu)))*nu), 0, 0, 0],
[float((E/((1+nu)*(1-2*nu)))*nu), float((E/((1+nu)*(1-2*nu)))*(1-nu)), float((E/((1+nu)*(1-2*nu)))*nu), 0, 0, 0],
[float((E/((1+nu)*(1-2*nu)))*nu), float((E/((1+nu)*(1-2*nu)))*nu), float((E/((1+nu)*(1-2*nu)))*(1-nu)), 0, 0, 0],
[ 0, 0, 0, float(E/(1+nu)), 0, 0],
[ 0, 0, 0, 0, float(E/(1+nu)), 0],
[ 0, 0, 0, 0, 0, float(E/(1+nu))]])
def epsilon(v):
return ufl.sym(ufl.grad(v))
def strain2voigt(e):
return ufl.as_vector([e[0,0],e[1,1],e[2,2],2*e[1,2], 2*e[0,2],2*e[0,1]])
def voigt2stress(s):
return ufl.as_tensor([
[s[0], s[5], s[4]],
[s[5], s[1], s[3]],
[s[4], s[3], s[2]]
])
def sigma(v,C):
return voigt2stress(ufl.dot(C, strain2voigt(epsilon(v))))
# Create mesh and define function space
V = fem.functionspace(domain, ("Lagrange", 1, (3,)))
bc_dofs = fem.locate_dofs_topological(V, 2, facets.find(101))
bcs = fem.dirichletbc(np.zeros((3,)), bc_dofs, V)
ds = ufl.Measure("ds", domain=domain,subdomain_data=facets)
dx = ufl.Measure("dx", domain=domain,subdomain_data=markers)
normal = ufl.FacetNormal(domain)
u = ufl.TrialFunction(V)
d = domain.geometry.dim
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0, 0)))
T = fem.Constant(domain,-1.0e6)
a = 0
start = time.time()
for i in range(n):
a += ufl.inner(sigma(u,C), epsilon(v))*dx(1000+i)
end = time.time()
print(end - start)
L = ufl.dot(f, v)*ufl.dx + ufl.dot(T*normal, v)*ds(102)
# Compute solution
#u = fem.Function(V)
start = time.time()
problem = LinearProblem(a, L, bcs=[bcs], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
end = time.time()
print(end - start)
start = time.time()
uh = problem.solve()
end = time.time()
print(end - start)
pyvista.start_xvfb()
# Create plotter and pyvista grid
p = pyvista.Plotter()
topology, cell_types, geometry = plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
# Attach vector values to grid and warp grid by vector
grid["u"] = uh.x.array.reshape((geometry.shape[0], 3))
actor_0 = p.add_mesh(grid, style="wireframe", color="k")
warped = grid.warp_by_vector("u", factor=1)
actor_1 = p.add_mesh(warped, show_edges=True)
p.show_axes()
if not pyvista.OFF_SCREEN:
p.show()
print(uh.vector.max())
return uh.x.array.max()
fe_case(4000)
In this code there is no C[i], so how am I to interpret this.
Please make an effort to describe your problem accurately.
@dokken thank you for being patient with me here is the corrected code:
import gmsh
import sys
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np
import matplotlib.pyplot as plt
import dolfinx
import time
def make_layers(n):
length = 1000;width = 1000; thickness= 0.5
gmsh.initialize()
gmsh.option.setNumber("General.Verbosity",0)
gmsh.model.add("plate")
for i in range(n):
gmsh.model.occ.addBox(0,0,i*thickness,length,width,thickness)
gmsh.model.occ.synchronize()
vols = gmsh.model.getEntities(3)
for i in range(n-1):
gmsh.model.occ.fragment([vols[i]],[vols[i+1]],removeObject=True,removeTool=True)
gmsh.model.occ.synchronize()
lines = gmsh.model.getEntities(1)
lengths = []
for i in lines:
lengths.append(gmsh.model.occ.getMass(1,i[1]))
surf = gmsh.model.getEntities(2)
vols = gmsh.model.getEntities(3)
for i in range(len(lines)):
if lengths[i] == length:
gmsh.model.mesh.setTransfiniteCurve(lines[i][1],2)
else:
gmsh.model.mesh.setTransfiniteCurve(lines[i][1],0)
for i in surf:
gmsh.model.mesh.setTransfiniteSurface(i[1])
gmsh.model.mesh.setRecombine(i[0],i[1])
for i in vols:
gmsh.model.mesh.setTransfiniteVolume(i[1])
f_surfaces = gmsh.model.getEntitiesInBoundingBox(-length/2,-width/2,-thickness/2,length/2,(n+1)*width,(n+1)*thickness,2)
l_surfaces = gmsh.model.getEntitiesInBoundingBox(length/2,-width/2,-thickness/2,2*length,(n+1)*width,(n+1)*thickness,2)
fixed_surfaces = []
load_surfaces = []
vol_count = 1000
for i in range(n):
fixed_surfaces.append(f_surfaces[i][1])
load_surfaces.append(l_surfaces[i][1])
gmsh.model.addPhysicalGroup(3,[i+1],vol_count)
gmsh.model.setPhysicalName(3,vol_count,"volume"+ str(i+1))
vol_count += 1
gmsh.model.addPhysicalGroup(2,fixed_surfaces,101)
gmsh.model.setPhysicalName(2,101,"fixed_surface")
#print(fixed_surfaces)
gmsh.model.addPhysicalGroup(2,load_surfaces,102)
gmsh.model.setPhysicalName(2,102,"force_surface")
#print(load_surfaces)
gmsh.model.mesh.generate(3)
gmsh.model.mesh.refine()
# if '-nopopup' not in sys.argv:
# gmsh.fltk.run()
domain, markers, facets = io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)
gmsh.finalize()
return domain,markers,facets
def fe_case(n):
domain,markers,facets = make_layers(n)
E = 1.1e9;nu = 0.3
stiffness_matrix = ufl.as_matrix([[float((E/((1+nu)*(1-2*nu)))*(1-nu)), float((E/((1+nu)*(1-2*nu)))*nu), float((E/((1+nu)*(1-2*nu)))*nu), 0, 0, 0],
[float((E/((1+nu)*(1-2*nu)))*nu), float((E/((1+nu)*(1-2*nu)))*(1-nu)), float((E/((1+nu)*(1-2*nu)))*nu), 0, 0, 0],
[float((E/((1+nu)*(1-2*nu)))*nu), float((E/((1+nu)*(1-2*nu)))*nu), float((E/((1+nu)*(1-2*nu)))*(1-nu)), 0, 0, 0],
[ 0, 0, 0, float(E/(1+nu)), 0, 0],
[ 0, 0, 0, 0, float(E/(1+nu)), 0],
[ 0, 0, 0, 0, 0, float(E/(1+nu))]])
C = stiffness_matrix*np.ones(n) # all stiffness matrices are same just for simplicity , ideally the values would be different for different subdomains
def epsilon(v):
return ufl.sym(ufl.grad(v))
def strain2voigt(e):
return ufl.as_vector([e[0,0],e[1,1],e[2,2],2*e[1,2], 2*e[0,2],2*e[0,1]])
def voigt2stress(s):
return ufl.as_tensor([
[s[0], s[5], s[4]],
[s[5], s[1], s[3]],
[s[4], s[3], s[2]]
])
def sigma(v,C):
return voigt2stress(ufl.dot(C, strain2voigt(epsilon(v))))
# Create mesh and define function space
V = fem.functionspace(domain, ("Lagrange", 1, (3,)))
bc_dofs = fem.locate_dofs_topological(V, 2, facets.find(101))
bcs = fem.dirichletbc(np.zeros((3,)), bc_dofs, V)
ds = ufl.Measure("ds", domain=domain,subdomain_data=facets)
dx = ufl.Measure("dx", domain=domain,subdomain_data=markers)
normal = ufl.FacetNormal(domain)
u = ufl.TrialFunction(V)
d = domain.geometry.dim
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0, 0)))
T = fem.Constant(domain,-1.0e6)
a = 0
start = time.time()
for i in range(n):
a += ufl.inner(sigma(u,C[i]), epsilon(v))*dx(1000+i)
end = time.time()
print(end - start)
L = ufl.dot(f, v)*ufl.dx + ufl.dot(T*normal, v)*ds(102)
# Compute solution
#u = fem.Function(V)
start = time.time()
problem = LinearProblem(a, L, bcs=[bcs], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
end = time.time()
print(end - start)
start = time.time()
uh = problem.solve()
end = time.time()
print(end - start)
pyvista.start_xvfb()
# Create plotter and pyvista grid
p = pyvista.Plotter()
topology, cell_types, geometry = plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
# Attach vector values to grid and warp grid by vector
grid["u"] = uh.x.array.reshape((geometry.shape[0], 3))
actor_0 = p.add_mesh(grid, style="wireframe", color="k")
warped = grid.warp_by_vector("u", factor=1)
actor_1 = p.add_mesh(warped, show_edges=True)
p.show_axes()
if not pyvista.OFF_SCREEN:
p.show()
print(uh.vector.max())
return uh.x.array.max()
fe_case(4000)
Make a single C where your np.ones(n) is replaced by a DG, 0 function containing a different constant per material.
Then you can write:
as
Q = dolfinx.fem.functionspace(mesh, ("DG",0))
D = dolfinx.fem.Function(Q)
# Assign i to each subdomain (100i)
for i in range(n):
D.interpolate(lambda x: np.full_like(x[0], i), cells=markers.find(1000+i))
C = stiffness_matrix * D
ids = tuple(1000+i for i in range(n))
a += ufl.inner(sigma(u,C), epsilon(v))*dx(ids)
this totally works.
small doubt: what exactly is x here and why x[0] ??
x is the input coordinate in the domain you are interpolating on, it is a vector (x[0], x[1], x[2]) where x[0] contains all x coordinates, x[1] all y etc.
I chose this to illustrate that if you want different values in the different regions you can do it in interpolate. In the current example, it would set the value in domain 1000+i to i. You can of course choose anything.
Got it . Thanks a lot!