Picking subdomains and assigning material properties for imported mesh

mesh
1

Hello, I’m trying to assign materials (and therefore picking subdomains) for imported mesh. I’m using an external tool to create mesh and then exporting it as med file format. I’m using meshio to convert mesh. I’m also exporting all meshes for all subdomains in order to further mark which element of whole domain (in file 1_full) belongs to which subdomain (files 1_air, 1_iron and 1_wire).

I’m doing it simply by assumption, that i might have domains of complex shape, so creating both geometry and mesh using external tool is preferable. The shape of those subdomains is complex enough, so that using any reasonable analytical function to mark boundary is impossible.

I’ve figured out, that in order to assign materials to elements in mesh 1_full, FEniCS is checking edges as well as cell centers. In order to find cell center I’ve created additional method that does exactly that.

This algorithm works, but it becomes computationally heavy and time-consuming for not that high number of elements (order of tens of thousands).

Is there any simpler way to choose subdomains and assign materials when dealing with imported mesh?

import meshio

print("Converting meshes...")
mesh = meshio.read('1_full.med')
meshio.write("1_full.xml", mesh)

mesh = meshio.read('1_air.med')
meshio.write("1_air.xml", mesh)

mesh = meshio.read('1_wire.med')
meshio.write("1_wire.xml", mesh)

mesh = meshio.read('1_iron.med')
meshio.write("1_iron.xml", mesh)

#############################################################################

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import mshr
import copy
from scipy import constants

print("Loading mesh...")

mesh = Mesh("1_full.xml")
mesh_iron = Mesh("1_iron.xml")
mesh_wire = Mesh("1_wire.xml")
mesh_air = Mesh("1_air.xml")

plot(mesh)
plt.show()

tol = 1e-14
def getMidPoint(pt1, pt2):
	return (pt1+pt2)*0.5

def getLine(pt1, pt2):
	A = -(pt1[1] - pt2[1])
	B = (pt1[0] - pt2[0])
	C = pt1[0]*A + pt1[1]*B
	return {'A': A, 'B': B, 'C': C}

def getMeshCellCenters(mesh):
	centers = []
	for cell in mesh.cells():
		pts = []
		for point in cell:
			pts.append(copy.copy(mesh.coordinates()[point]))
			
		line1 = getLine(pts[0], getMidPoint(pts[1], pts[2]))
		line2 = getLine(pts[1], getMidPoint(pts[0], pts[2]))
		
		W = line1['A']*line2['B'] - line2['A']*line1['B']
		Wx = line1['C']*line2['B'] - line2['C']*line1['B']
		Wy = line1['A']*line2['C'] - line2['A']*line1['C']
		centers.append([Wx/W, Wy/W])
		
	return np.array(centers)

class MeshDomain(SubDomain):
	def __init__(self):
		SubDomain.__init__(self)
		self.nodes = []
		self.cell_centers = []
		
	def physicalSurface(self, mesh):
		self.nodes = mesh.coordinates()
		self.cell_centers = getMeshCellCenters(mesh)

	def inside(self, point, on_boundry):
		for pt in self.nodes:
			if near(point[0], pt[0], tol) & near(point[1], pt[1], tol):
				return True
		for pt in self.cell_centers:
			if near(point[0], pt[0], tol) & near(point[1], pt[1], tol):
				return True
		return False
			
print("Assigning materials...")

materials = MeshFunction('size_t', mesh, mesh.topology().dim())

subdom_air = MeshDomain()
subdom_iron = MeshDomain()
subdom_wire = MeshDomain()

subdom_air.physicalSurface(mesh_air)
subdom_iron.physicalSurface(mesh_iron)
subdom_wire.physicalSurface(mesh_wire)

subdom_air.mark(materials, 0)
subdom_iron.mark(materials, 1)
subdom_wire.mark(materials, 2)

print("Saving materials to file...")
File('materials.xml') << materials
2

Hi, consider the following

from scipy.spatial import cKDTree
from dolfin import *
import numpy as np


mesh = UnitSquareMesh(256, 256)

tdim = mesh.topology().dim()
# A piece of mesh where we have the material properties
f = MeshFunction('size_t', mesh, tdim, 0)
CompiledSubDomain('x[0] < 0.5 + DOLFIN_EPS').mark(f, 1)
submesh = SubMesh(mesh, f, 1)

# Submesh has a mapping back to mesh cells but we want to compute it
# in case where Mesh is not a SubMesh instance
parent = np.array([c.midpoint().array() for c in cells(mesh)])
child = np.array([c.midpoint().array() for c in cells(submesh)])

tree = cKDTree(parent)
_, child2parent = tree.query(child, k=1)

# It is correct
child2parent0 = submesh.data().array('parent_cell_indices', tdim)

assert np.linalg.norm(child2parent - child2parent0) < 1E-15

fsub = MeshFunction('size_t', submesh, tdim, 1)
# How to use it to assign to CellFunction on
f = MeshFunction('size_t', mesh, tdim, 0)
f.array()[child2parent] = fsub.array()

a = Constant(1)*dx(domain=mesh, subdomain_data=f, subdomain_id=1)
assert near(assemble(a), 0.5)
1 Like