Assigning different physical groups to different color in gmsh and mark them as different material in fenics

General
29

@dokken Yeah sure, I will share the code here but I am not getting the results I wanted (please see the image attached) so kindly help that in which step I making mistake, here’s the code:

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

V = FunctionSpace(mesh, "DG", 0)
u = Function(V)

values = np.arange(4383*1) 
local_values = np.zeros_like(u.vector().get_local())
local_values_material = np.zeros_like(u.vector().get_local())

for cell in cells(mesh):
    midpoint = cell.midpoint().array()
    i = int(midpoint[0])
    j = int(midpoint[1])
    k = int(midpoint[2])
    ## partha idea
    local_values_material[cell.index()] = material_array_final[cell.index()]
    materials[cell] = int(local_values_material[cell.index()])
    print(midpoint, i, j , k, "Material:", materials[cell], "Cell index: ", cell.index())
    
u.vector().set_local(local_values)
dolfin.XDMFFile(dolfin.MPI.comm_world, "mesh1.xdmf").write_checkpoint(u,"u",0)

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] == 1:
            values[0] = self.k_0
        else:
            values[0] = self.k_1
            
E1 = 210e3
nu1 = 0.3
E2 = 21e3
nu2 = 0.25

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

File('e.pvd') << project(E, FunctionSpace(mesh, 'DG', 0))
File('materials.pvd') << materials

Screenshot 2023-06-04 161859
Screenshot 2023-06-04 1618591913×1077 451 KB

I wanted the material properties as shown in gmsh but I got the right one as shown in paraview. I’m not able to figure it out that when I have assigned each cell with a cell tag of either 0 or 1 and calling the same in user-defined function then where am I doing mistake.
Second Question: I tried making another xdmf file as @dokken suggested for facets by referring to this post, but was not able to make it.
I tried:

import meshio

msh = meshio.read("/mnt/d/Research Projects/FEM for microstructures/sample1.msh")

meshio.write("mesh.xdmf",
    meshio.Mesh(points=msh.points,
        cells={"triangle": msh.cells_dict["triangle"]},
    )
)

meshio.write("mf.xdmf", meshio.Mesh(points=msh.points[:,:2], cells={"line": msh.cells_dict["line"]}))

from dolfin import * 

mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

mvc = MeshValueCollection("size_t", mesh, 1) 
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

plot(mesh)

got the error:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/tmp/ipykernel_79/2720985724.py in <module>
     19 mvc = MeshValueCollection("size_t", mesh, 1)
     20 with XDMFFile("mf.xdmf") as infile:
---> 21     infile.read(mvc, "name_to_read")
     22 mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
     23 

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to recognise cell type.
*** Reason:  Unknown value "".
*** Where:   This error was encountered inside XDMFFile.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  ubuntu
*** -------------------------------------------------------------------------

Any input will be appreciated.
Thank You