chfem Python API Tutorial

# install chfem
!pip install git+https://gitlab.com/cortezpedro/chfem_gpu.git@dev -q

# To enable c stdout printing in colab notebook
!pip install wurlitzer -q
import wurlitzer
wurlitzer.Wurlitzer.flush_interval=0.001
%load_ext wurlitzer

# necessary for pyvista plots
!pip install pyvista -q
!pip install piglet -q
!pip install pyvirtualdisplay -q
!apt-get -qq install xvfb
from pyvirtualdisplay import Display
display = Display(visible=0, size=(600, 400))
display.start()
def pv_plot(array, cmap='jet', opacity=None):
  grid = pv.ImageData()
  grid.dimensions = ( array.shape[2], array.shape[1], array.shape[0] )
  grid.point_data["scalars"] = np.reshape( array, (array.size) )
  plotter = pv.Plotter(notebook=True)
  if opacity is None:
    plotter.add_volume(grid, scalars="scalars", cmap=cmap)
  else:
    plotter.add_volume(grid, scalars="scalars", cmap=cmap, opacity=opacity)
  plotter.show()

# other imports
from matplotlib import pyplot as plt
import pyvista as pv
import numpy as np
import chfem

def create_cylinder(shape, radius, position, direction='z', clr=255, array=None):
  if array is None:
    array = np.zeros(shape,dtype='uint8')
  if direction.lower() == 'x':
    array = array.transpose(2,0,1) # x z y
    return create_cylinder(array.shape,radius,position,'z',clr,array).transpose(1,2,0)
  elif direction.lower() == 'y':
    array = array.transpose(1,2,0) # y x z
    return create_cylinder(array.shape,radius,position,'z',clr,array).transpose(2,0,1)
  if not (direction.lower() == 'z'):
    raise ValueError(f'unexpected direction: {direction}')
  x,y = position
  dy2_arr = np.power(np.arange(shape[1])-y+0.5,2.)
  for xp in range( max(0,int(np.floor(x-radius))), min(int(np.ceil(x+radius)+1),shape[2]) ):
    dist = np.sqrt( (xp+0.5-x)**2 + dy2_arr )
    array[ 0, dist<=radius , xp ] = clr
  array[1:,array[0]==clr]=clr
  return array

shape = (300,150,200) # (z,y,x) == (slices,rows,cols)
array = create_cylinder(shape,20,(100,75),'z')
array = create_cylinder(array.shape,20,(75,150),'x',array=array)
array = create_cylinder(array.shape,20,(150,100),'y',array=array)

pv_plot(array)

keff = chfem.compute_property('conductivity', array, mat_props={255: 1, 0: 0.1}, direction='x', output_fields="cubes")

temperature = chfem.import_scalar_field_from_chfem("cubes_temperature_0.bin", array.shape)
pv_plot(temperature)

Ceff = chfem.compute_property('elasticity', array, mat_props={255: (200, 0.2), 0: (100, 0.3)}, direction='x', output_fields="cubes")

displacement = chfem.import_vector_field_from_chfem("cubes_displacement_0.bin", array.shape)
pv_plot(displacement[0]) # x direction

def gaussian_opacity(img, sigma_multiplier=1., n=101):
  mu = np.mean(img)
  sigma = np.std(img)
  sigma *= sigma_multiplier
  return 1. - np.exp( -( np.linspace(np.min(img),np.max(img),n) -mu  )**2 / (2*sigma*sigma) )

# making the phase in 255 be the pores
if 255 in array:
  array[array==0]=1
  array[array==255]=0

Keff = chfem.compute_property('permeability', array, direction='x', output_fields="cubes")

pressure = chfem.import_scalar_field_from_chfem("cubes_pressure_0.bin", array.shape)
pv_plot(pressure, opacity=gaussian_opacity(pressure))