How to set up a function space consisting of multiple domains from an imported gmsh mesh for an elasticity analysis

[prusso630]

I have an imported gmsh consisting of two domains (name: slab,tag: 11) and (name: rebar, tag: 12).

from contextlib import ExitStack

import numpy as np


from dolfinx import la
from dolfinx.fem import (Expression, Function, FunctionSpace,
                         VectorFunctionSpace, dirichletbc, form,
                         locate_dofs_topological)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
                               set_bc)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_box,
                          locate_entities_boundary)
from ufl import dx, grad, inner

from mpi4py import MPI
from petsc4py import PETSc

dtype = PETSc.ScalarType


from mpi4py import MPI

import meshio 


## original mesh file
msh = meshio.read("mesh3D.msh")

## physical surface & volume data
for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]
    elif key == "tetra":
        tetra_data = msh.cell_data_dict["gmsh:physical"][key]

## cell data
tetra_cells = np.array([None])
triangle_cells = np.array([None])
for cell in msh.cells:
    if cell.type == "tetra":
        if tetra_cells.all() == None:
            tetra_cells = cell.data
        else:
            tetra_cells = np.concatenate((tetra_cells,cell.data))
    elif cell.type == "triangle":
        if triangle_cells.all() == None:
            triangle_cells = cell.data
        else:
            triangle_cells = np.concatenate((triangle_cells,cell.data))



print('done')

The lengths for tetra_cells and tetra_data are the same. tetra_data seems to contain the group tags 11 and 12 for each cell.

There will be a modulus of elasticity for group 11 and 12 for concrete and steel respectively to account for. What I would like to start to do at first is set up a function space system for the two domains that are imported. There will be a boundary condition at nodes near Z=-72 and Z=72 where displacement should be zero for a fixed support at each end. I guess to set up an analysis is maybe beyond the scope of this question so far I just would like to see how to start to set that up…

Anyone can help me please how to set up for the two imported domains a function space taking into account the two modulus E’s and the boundary conditions?

P.S. My current version for dolfinx is V0.6.0 installed with spack. I tried to use the import function from:

Subdomain for different materials based on physical volume

@dokken Thank you for the suggestion. I have solved my issue, but for other people who might be looking for a tip (that are not trivial to beginners), I will describe my issue and the solution below: After hours of investigation, I have noticed there were three main issues I’ve had: 1: xdmf output My model is consisted of two objects, and each of them are stored as different elements in msh.cells array. This means there are two separated cell data sections whose type is ‘tetra’ (and same goes…

however I could not yet get it moved onto dolfinx v0.6.0 there were some import challenges to get through where imports have changed. Anyone can help me to find where those are in dolfinx v0.6…0?

[dokken]

I would suggest reading: Defining subdomains for different materials — FEniCSx tutorial for material parameter specification, and Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial for boundary conditions

[prusso630]

import gmsh
import numpy as np
import pyvista
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities,  locate_entities_boundary 
from dolfinx.plot import create_vtk_mesh

from ufl import (SpatialCoordinate, TestFunction, TrialFunction,
                 dx, grad, inner)

from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from dolfinx.io import gmshio


mesh, cell_markers, facet_markers = gmshio.read_from_msh("mesh3D.msh", MPI.COMM_WORLD, gdim=3) 

Q = FunctionSpace(mesh, ("DG", 1))



E = Function(Q)
#cell_markers0 = locate_entities_boundary(mesh, mesh.topology.dim, 11)
#cell_markers1 = locate_entities_boundary(mesh, mesh.topology.dim, 12)


print('done')

I see that cell_markers is an array of the integer values 11 and 12. Following “Defining subdomain for different materials” I was able to get as far as the sample above. As opposed to locating entities by co-ordinate value I have to locate them by way of cell_marker’s array. I found locate_entities_boundary nearrby and was thinking maybe it can be used it.

There is a function defenition for it:

def locate_entities_boundary(mesh: Mesh, dim: int, marker: typing.Callable) -> np.ndarray:

The below code does not work I think because locate_entities_boundary expects a typing.Callable…

    cell_markers0 = locate_entities_boundary(mesh, mesh.topology.dim, 11)

*I was trying to start to build a map in this case of mu, lambda where each respectively is an array of petsc type scalar which will vary by the different elasticity that belong to domain 11 and domain 12 as listed in cell_markers…

Each entry in mu’s and lambda’s array will be unique by the elasticity of the domain.*

How can the entities on domains 11 and 12 listed in cell_markers be located?
(As in similar fashion to:)

def Omega_0(x):
    return x[1] <= 0.5

def Omega_1(x):
    return x[1] >= 0.5

kappa = Function(Q)
cells_0 = locate_entities(mesh, mesh.topology.dim, Omega_0)
cells_1 = locate_entities(mesh, mesh.topology.dim, Omega_1)

P.S.::::
I have so far:

from dolfinx.io import gmshio


mesh, cell_markers, facet_markers = gmshio.read_from_msh("mesh3D.msh", MPI.COMM_WORLD, gdim=3) 

Q = FunctionSpace(mesh, ("DG", 1))



E = Function(Q)
mu_indices = np.where(cell_markers.values == 11)[0]
lambda_indices = np.where(cell_markers.values == 12)[0]
E.x.array[mu_indices] = np.full_like(mu_indices, 200, dtype=ScalarType)
E.x.array[lambda_indices] = np.full_like(lambda_indices, 500, dtype=ScalarType)

So far the len of E.x.array is showing too large if I am not mistaken. There should be:

len(cell_markers.values)
19690
len(E.x.array)
78760 <<--- looks too big to me there 

The data type of array mu_indices, lambda_indices is not yet correct. It was the best I could do so far…
I see that there are a meshtags and meshtags_from_entities but I am not sure yet about the usage of those if those can be of help in this situation…

[prusso630]

I also tried from:

Subdomain for different materials based on physical volume

Hi, I hope there is an existing answer: I have a complicated geometry with multiple materials combined, and I wish to use constant values that are unique to each material (e.g. Kappa in https://fenicsproject.org/pub/tutorial/sphinx1/._ftut1005.html). These materials are defined as Physical Volume in the mesh, so I wish to extract this data and apply that information to define these constants (k=1 for Physical Volume 1, k=1.25 for Physical Volume 2, etc.) without definind the boundary with mat…

import gmsh
import numpy as np
import pyvista
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square, locate_entities,  locate_entities_boundary
from dolfinx.plot import create_vtk_mesh
from dolfinx.mesh import *
from ufl import (SpatialCoordinate, TestFunction, TrialFunction,
                 dx, grad, inner)
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from dolfinx.io import gmshio


mesh, materials, facet_markers = gmshio.read_from_msh("mesh3D.msh", MPI.COMM_WORLD, gdim=3) 

class EMod(UserExpression):
    def __init__(self,materials,volume_list,rho_list,**kwargs):
        super().__init__(**kwargs)
        self.materials = materials
        self.volume_list = volume_list
        self.rho_list = rho_list
    def value_shape(self):
        return ()
    def eval_cell(self,values,x,cell):
        values[0] = 0
        for i in range(len(self.volume_list)):
            if self.materials[cell.index] == self.volume_list[i]:
                values[0] = self.rho_list[i]   

volume_list = [11,12]                          # list of physical volume ID's
rho_list    = [100,200] # list of pre-defined values
rho = EMod(materials,volume_list,rho_list,degree=1)


#E.x.array[mu_indices] = np.full_like(mu_indices, 200, dtype=ScalarType)
#E.x.array[lambda_indices] = np.full_like(lambda_indices, 500, dtype=ScalarType)
print('done')

The difficulty here is that UserExpression is not being found. I looked around about UserExpression and found it to be deprected. Is it possible to switch over directly to Expression as opposed to UserExpression which was from fenics?

[dokken]

@prusso630

Q = FunctionSpace(mesh, ("DG", 1))



E = Function(Q)
mu_indices = np.where(cell_markers.values == 11)[0]
lambda_indices = np.where(cell_markers.values == 12)[0]
E.x.array[mu_indices] = np.full_like(mu_indices, 200, dtype=ScalarType)
E.x.array[lambda_indices] = np.full_like(lambda_indices, 500, dtype=ScalarType)

You are assuming a 1 to 1 relationship between the dofs of a DG-1 space and the cell indices. This is not correct, as there is only a 1-1 relationship between a DG-0 space and the cell indices. So i would suggest using a DG-0 space.

See also part two of the tutorial: Defining subdomains for different materials — FEniCSx tutorial
which creates materials based on cell data from the mesh markers.