Heat Transfer using custom image

mesh
1

Hello everyone in this community,

I have been trying to solve a heat transfer problem using my image which I project on the mesh.
I followed the method explained in : Manipulating images
The image array contains a property and when I plot it I get as follows. The one on the left is is correct and what I want:

combined
combined1279Ɨ545 88.4 KB

This image is based on the following method:
l=1

mesh = dl.RectangleMesh(dl.Point(0.,0.),dl.Point(1.*l,1.*l), nx, ny)

ph = np.zeros((nx,ny))
def pixel():
for i in range(len(markers[:,0])):
for j in range(len(markers[0,:])):
for k in dic:
if markers[i][j] == k:
ph[i][j] = 1/(dic[k])
break
else:
ph[i][j] = 0.048

return ph

pixel()

class FE_image(UserExpression):
def eval_cell(self, value, x, ufc_cell):
p = Cell(mesh, ufc_cell.index).midpoint()
i, j = int(p[0](nx)), int(p[1](ny))
value[:] = ph[-(j+1), i]

def value_shape(self):
    return ()

y = FE_image()
print(y)

V = FunctionSpace(mesh, ā€˜CGā€™, 1)
u0 = project(y, V)
k_func = Function(V)
k_func.assign(u0)

Here ph is the property array and then using the method mentioned in the previous link to assign them to the mesh and later on project it to the mesh space V. k_func is the space with the values which I plot as shown in the previous image.

Now, for my simulation the domain size of 1x1 is pretty big and so I decrease the value of l from 1 to 7e-2 but this resulted in a weird plot of k_func and captured only a part of it which you can see in the right image

I thought this method would first create a mesh with a domain size and number of elements and then assign the array to the mesh but it is not doing that.
I went through other threads but it was pretty different.
Can someone please guide with this or have any other idea of using an image in fenics?

Thank You!

2

Image of the code for easy understanding:

code
code632Ɨ679 52.9 KB

3

First of all, please use 3x` encapsulation to make your code copyable, Ref

Please also make sure that the complete code is added (including definitions of your image).

And finally if you scale your mesh with l, then this has to be taken into account when computing i and j

4

Hello @dokken
Thank you for replying. I am sorry and I am providing the code using the preformatted text.

#Import Modules
from __future__ import print_function
from fenics import *
import ufl
import numpy as np
import matplotlib.pyplot as plt
import cv2
import pandas as pd
from PIL import Image
from scipy import ndimage
import dolfin as dl

#Image Input
img1 = cv2.imread("sample_fea_mod.png")
img = cv2.cvtColor(img1,cv2.COLOR_BGR2GRAY)


#Thresholding and labeling
ret,thresh = cv2.threshold(img,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
ret3, markers = cv2.connectedComponents(thresh)
(nx, ny) = markers.shape

k_s= 0.8
k_bf = 0.5  
beta = 1     
l1 = 100
rho_s= 2650
c_s= 700


props = measure.regionprops_table(markers, intensity_image=img,
                                  properties = ['label',
                                                'centroid',
                                               'area',
                                                'equivalent_diameter', 
                                                'perimeter'])   


df = pd.DataFrame(props)

df = df.drop(0)
print(df)
df.to_csv('file_name.csv')
label_list = list(df["label"])
dia_list1 = list(df["equivalent_diameter"])
dia_list = [item * pixels_to_um for item in dia_list1]
dic = dict(zip(label_list, dia_list))


#Mesh
l=1

mesh = dl.RectangleMesh(dl.Point(0.,0.),dl.Point(1.*l,1.*l), nx, ny)


#Property Array
ph = np.zeros((nx,ny))
def pixel():
    for i in range(len(markers[:,0])):
        for j in range(len(markers[0,:])):
            for k in dic:
                if markers[i][j] == k:
                    ph[i][j] = 1/(dic[k])
                    break
                else:
                    ph[i][j] = 0.048
                
    return ph

pixel()


#Mapping array
class FE_image(UserExpression):
    def eval_cell(self, value, x, ufc_cell):
        p = Cell(mesh, ufc_cell.index).midpoint() 
        i, j = int(p[0]*(nx)), int(p[1]*(ny))
        value[:] = ph[-(j+1), i]

    def value_shape(self):
        return ()

y = FE_image()
print(y)

V = FunctionSpace(mesh, 'CG', 1)
u0 = project(y, V)
k_func = Function(V)
k_func.assign(u0)

kf = File("cond_new.pvd")
kf<<k_func

So, k_func has my conductivity values projected in the function space which I am saving. Now l = 1 gives me a unit square mesh and my output image is correct. Decrease l to 7e-2 gives the wrong output as shown in the right figure.

5

As I said in my previous post,

you have to adjust i and j as the length/height of the domain is no longer 1.

should be

    i, j = int(p[0]/l*(nx)), int(p[1]/l*(ny))`
1 Like
6

This works.
Thank you so much!