diff --git a/codes/.vscode/launch.json b/codes/.vscode/launch.json new file mode 100644 index 0000000..17e15f2 --- /dev/null +++ b/codes/.vscode/launch.json @@ -0,0 +1,15 @@ +{ + // Use IntelliSense to learn about possible attributes. + // Hover to view descriptions of existing attributes. + // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 + "version": "0.2.0", + "configurations": [ + { + "name": "Python: Current File", + "type": "python", + "request": "launch", + "program": "${file}", + "console": "integratedTerminal" + } + ] +} \ No newline at end of file diff --git a/codes/.vscode/tasks.json b/codes/.vscode/tasks.json new file mode 100644 index 0000000..6c28385 --- /dev/null +++ b/codes/.vscode/tasks.json @@ -0,0 +1,12 @@ +{ + // See https://go.microsoft.com/fwlink/?LinkId=733558 + // for the documentation about the tasks.json format + "version": "2.0.0", + "tasks": [ + { + "label": "echo", + "type": "shell", + "command": "echo Hello" + } + ] +} \ No newline at end of file diff --git a/codes/CS.py b/codes/CS.py new file mode 100644 index 0000000..1d0d2bd --- /dev/null +++ b/codes/CS.py @@ -0,0 +1,705 @@ +import numpy as np +from numpy import linalg as LA +import sys +from mpi4py import MPI +comm = MPI.COMM_WORLD +size = comm.Get_size() +rank = comm.Get_rank() + + +# COMPRESSED SENSING: LINEAR BREGMAN METHOD +# Translated and adapted into python from tinycs +# +# *tinycs* is a minimal compressed sensing (CS) toolkit designed +# to allow MR imaging scientists to design undersampled +# acquisitions and reconstruct the resulting data with CS without +# needing to be a CS expert. +# +# The Cartesian reconstruction is based on the split Bregman +# code written by Tom Goldstein, originally available here: +# + + +def pdf(k,kw,klo,q): + + p = (np.abs(k)/kw)**(-q) + p[np.where(k==0)] = 0 + p[np.where(np.abs(k)<=kw)] = 1 + p[np.where(k=norm: + break + + P = np.fft.fftshift(P) + + + #if np.mod(n,2)!=0: + # P = np.concatenate(([1],P),axis=None) + + + + return P + +def mask_pdf_2d(dims,norm,q,pf): + + nz = dims[1] + ny = dims[0] + yc = round(ny/2) + zc = round(nz/2) + rmax = np.sqrt((ny-yc)**2 + (nz-zc)**2) + [Z,Y] = np.meshgrid(np.arange(0,nz),np.arange(0,ny)) + RR = np.sqrt( (Y-yc)**2 + (Z-zc)**2) + Z = np.abs(Z - nz/2 - 0.5) + Y = np.abs(Y - ny/2 - 0.5) + + + for rw in range(1,int(rmax)+1): + P = np.ones([ny,nz])/pf + C = np.logical_and( Z <= rw , Y <= rw) + W = np.logical_or( Z > rw , Y > rw) + P[W] = (RR[W]/rw)**(-q) + if np.sum(P) >= norm: + break + + + + return [P,C] + +def GeneratePattern(dim,R): + + # 3D CASE + if np.size(dim)==3: + + nro = dim[0] + npe = dim[1] + nacq = round(npe/R) + q = 1 + pf = 1 + P = mask_pdf_1d(npe, nacq, q, pf) + while True: + M = np.random.rand(npe) + M = 1*(M<=P) + if np.sum(M)==nacq: + break + # remove partial Fourier plane and compensate sampling density + M = M!=0 + M = np.tile(M,[nro,1]); + #M = M.T + + # 4D CASE + if np.size(dim)==4: + + nro = dim[0] + npe1 = dim[1] + npe2 = dim[2] + nacq = round(npe1*npe2/R) + q = 1 + pf = 1 + [P,C] = mask_pdf_2d([npe1,npe2], nacq, q, pf) + RR = np.random.rand(npe1,npe2) + M = (RR <= P) + nchosen = np.sum(M) + if nchosen > nacq: # Correct for inexact number chosen + #outerOn = np.logical_and( M , P!=1 ) + outerOn = np.where((M)*(P!=1)) + numToFlip = nchosen-nacq + idxs = np.random.permutation(outerOn[0].size) + idxx = outerOn[0][idxs[0:numToFlip]] + idxy = outerOn[1][idxs[0:numToFlip]] + M[idxx,idxy] = False + + elif nchosen < nacq: + outerOff = np.where(~M) + idxs = np.random.permutation(outerOff[0].size) + numToFlip = nacq - nchosen + idxx = outerOff[0][idxs[0:numToFlip]] + idxy = outerOff[1][idxs[0:numToFlip]] + M[idxx,idxy] = True + + + + + M = np.rollaxis(np.tile(np.rollaxis(M,1),[nro,1,1]),2) + M = np.fft.ifftshift(M) + M = M.transpose((1,0,2)) + + + return M + +def get_norm_factor(MASK,uu): + UM = MASK==1 + return UM.shape[0]/LA.norm(uu) + +def Dxyzt(X): + + if np.ndim(X)==3: + dd0 = X[:,:,0] + dd1 = X[:,:,1] + DA = dd0 - np.vstack((dd0[1::,:],dd0[0,:])) + DB = dd1 - np.hstack((dd1[:,1::],dd1[:,0:1])) + + return DA + DB + + if np.ndim(X)==4: + dd0 = X[:,:,:,0] + dd1 = X[:,:,:,1] + dd2 = X[:,:,:,2] + + DA = dd0 - np.vstack((dd0[1::,:,:],dd0[0,:,:][np.newaxis,:,:])) + DB = dd1 - np.hstack((dd1[:,1::,:],dd1[:,0,:][:,np.newaxis,:])) + DC = dd2 - np.dstack((dd2[:,:,1::],dd2[:,:,0][:,:,np.newaxis])) + + return DA + DB + DC + +def Dxyz(u): + + if np.ndim(u)==2: + + dx = u[:,:]- np.vstack((u[-1,:],u[0:-1,:])) + dy = u[:,:]- np.hstack((u[:,-1:],u[:,0:-1])) + D = np.zeros([dx.shape[0],dx.shape[1],2],dtype=complex) + D[:,:,0] = dx + D[:,:,1] = dy + + return D + + if np.ndim(u)==3: + + dx = u[:,:,:]- np.vstack((u[-1,:,:][np.newaxis,:,:],u[0:-1,:,:])) + dy = u[:,:,:]- np.hstack((u[:,-1,:][:,np.newaxis,:],u[:,0:-1,:])) + dz = u[:,:,:]- np.dstack((u[:,:,-1][:,:,np.newaxis],u[:,:,0:-1])) + + D = np.zeros([dx.shape[0],dx.shape[1],dx.shape[2],3],dtype=complex) + + D[:,:,:,0] = dx + D[:,:,:,1] = dy + D[:,:,:,2] = dz + + return D + +def shrink(X,pgam): + p = 1 + s = np.abs(X) + tt = pgam/(s)**(1-p) + # t = pgam/np.sqrt(s) + ss = s-tt + ss = ss*(ss>0) + s = s + 1*(s0.001) + + if scantype=='0G': + PHASE0[:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,k])>0.001) + 10*varPHASE0 + PHASE1[:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,k])>0.001) + 10*varPHASE0 + np.pi*Sq[:,:,k]/VENC + + if scantype=='-G+G': + PHASE0[:,:,k] = gamma*B0*TE*np.ones([row,col]) + 10*varPHASE0 - np.pi*Sq[:,:,k]/VENC + PHASE1[:,:,k] = gamma*B0*TE*np.ones([row,col]) + 10*varPHASE0 + np.pi*Sq[:,:,k]/VENC + + RHO0[:,:,k] = modulus*np.cos(PHASE0[:,:,k]) + Drho + 1j*modulus*np.sin(PHASE0[:,:,k]) + 1j*Drho2 + RHO1[:,:,k] = modulus*np.cos(PHASE1[:,:,k]) + Drho + 1j*modulus*np.sin(PHASE1[:,:,k]) + 1j*Drho2 + + + if np.ndim(Sq)==4: + [row, col, dep, numt2] = Sq.shape + [X, Y, Z] = np.meshgrid(np.linspace(0, col, col), np.linspace( + 0, row, row), np.linspace(0, dep, dep)) + + for k in range(numt2): + + if noise: + Drho = np.random.normal(0, 0.2, [row, col, dep]) + Drho2 = np.random.normal(0, 0.2, [row, col, dep]) + else: + Drho = np.zeros([row, col, dep]) + Drho2 = np.zeros([row, col, dep]) + + varPHASE0 = np.random.randint(-10, 11, size=(row, col, dep)) * \ + np.pi/180*(np.abs(Sq[:, :, :, k]) < 0.001) + modulus = 0.5 + 0.5*(np.abs(Sq[:, :, :, k]) > 0.001) + + if scantype == '0G': + PHASE0[:, :, :, k] = (gamma*B0*TE+0.01*X) * \ + (np.abs(Sq[:, :, :, k]) > 0.001) + 10*varPHASE0 + PHASE1[:, :, :, k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:, :, :, k]) + > 0.001) + 10*varPHASE0 + np.pi*Sq[:, :, :, k]/VENC + + if scantype == '-G+G': + PHASE0[:, :, :, k] = gamma*B0*TE * \ + np.ones([row, col, dep]) + varPHASE0 - np.pi*Sq[:, :, :, k]/VENC + PHASE1[:, :, :, k] = gamma*B0*TE * \ + np.ones([row, col, dep]) + varPHASE0 + np.pi*Sq[:, :, :, k]/VENC + + RHO0[:, :, :, k] = modulus*np.cos(PHASE0[:, :, :, k]) + \ + Drho + 1j*modulus*np.sin(PHASE0[:, :, :, k]) + 1j*Drho2 + RHO1[:, :, :, k] = modulus*np.cos(PHASE1[:, :, :, k]) + \ + Drho + 1j*modulus*np.sin(PHASE1[:, :, :, k]) + 1j*Drho2 + + + + return [RHO0,RHO1] + +def undersampling(Sqx,Sqy,Sqz,options,savepath): + + R = options['cs']['R'] + + for r in R: + + if rank==0: + print('Using Acceleration Factor R = ' + str(r)) + print('Component x of M0') + + [M0,M1] = GenerateMagnetization(Sqx,options['cs']['VENC'],options['cs']['noise']) + + + print('\n Component x of M0') + M0_cs = CSMETHOD(M0,r) + print('\n Component x of M1') + M1_cs = CSMETHOD(M1,r) + + Sqx_cs = phase_contrast(M1_cs,M0_cs,options['cs']['VENC']) + del M0,M1 + del M0_cs, M1_cs + + [M0,M1] = GenerateMagnetization(Sqy,options['cs']['VENC'],options['cs']['noise']) + + + print('\n Component y of M0') + M0_cs = CSMETHOD(M0,r) + print('\n Component y of M1') + M1_cs = CSMETHOD(M1,r) + + + Sqy_cs = phase_contrast(M1_cs,M0_cs,options['cs']['VENC']) + + del M0,M1 + del M0_cs, M1_cs + + [M0,M1] = GenerateMagnetization(Sqz,options['cs']['VENC'],options['cs']['noise']) + + if rank==0: + print('\n Component z of M0') + M0_cs = CSMETHOD(M0,r) + if rank==0: + print('\n Component z of M1') + M1_cs = CSMETHOD(M1,r) + if rank==0: + print(' ') + + Sqz_cs = phase_contrast(M1_cs,M0_cs,options['cs']['VENC']) + + + if rank==0: + print('saving the sequences in ' + savepath) + seqname = options['cs']['name'] +'_R' + str(r) + '.npz' + print('sequence name: ' + seqname) + np.savez_compressed( savepath + seqname, x=Sqx_cs, y=Sqy_cs,z=Sqz_cs) + + del Sqx_cs,Sqy_cs,Sqz_cs + +def undersampling_peakpv(Sqx,Sqy,Sqz,options,R): + + Sqx_cs = {} + Sqy_cs = {} + Sqz_cs = {} + [Mx0,Mx1] = GenerateMagnetization(Sqx,options['cs']['VENC'],options['cs']['noise'],scantype='0G') + [My0,My1] = GenerateMagnetization(Sqy,options['cs']['VENC'],options['cs']['noise'],scantype='0G') + [Mz0,Mz1] = GenerateMagnetization(Sqz,options['cs']['VENC'],options['cs']['noise'],scantype='0G') + + + Mx0_cs = CSMETHOD(Mx0,R) + Mx1_cs = CSMETHOD(Mx1,R) + My0_cs = CSMETHOD(My0,R) + My1_cs = CSMETHOD(My1,R) + Mz0_cs = CSMETHOD(Mz0,R) + Mz1_cs = CSMETHOD(Mz1,R) + + Sqx_cs = phase_contrast(Mx1_cs,Mx0_cs,options['cs']['VENC'],scantype='0G') + Sqy_cs = phase_contrast(My1_cs,My0_cs,options['cs']['VENC'],scantype='0G') + Sqz_cs = phase_contrast(Mz1_cs,Mz0_cs,options['cs']['VENC'],scantype='0G') + + + + return [Sqx_cs,Sqy_cs,Sqz_cs] + +def undersampling_short(Mx,My,Mz,options): + + R = options['cs']['R'] + savepath = options['cs']['savepath'] + + + R_SENSE = 1 + if 'R_SENSE' in options['cs']: + R_SENSE = options['cs']['R_SENSE'][0] + + for r in R: + if rank==0: + print('Using Acceleration Factor R = ' + str(r)) + + + if R_SENSE==2: + [MxS0_cs,MxS1_cs] = CSMETHOD_SENSE(Mx,r,2) + [MyS0_cs,MyS1_cs] = CSMETHOD_SENSE(My,r,2) + [MzS0_cs,MzS1_cs] = CSMETHOD_SENSE(Mz,r,2) + if rank==0: + print('saving the sequences in ' + savepath) + seqname_s0 = options['cs']['name'] +'S0_R' + str(r) + '.npz' + seqname_s1 = options['cs']['name'] +'S1_R' + str(r) + '.npz' + print('sequence name: ' + seqname_s0) + np.savez_compressed( savepath + seqname_s0, x=MxS0_cs, y=MyS0_cs,z=MzS0_cs) + print('sequence name: ' + seqname_s1) + np.savez_compressed( savepath + seqname_s1, x=MxS1_cs, y=MyS1_cs,z=MzS1_cs) + del MxS0_cs, MyS0_cs, MzS0_cs + del MxS1_cs, MyS1_cs, MzS1_cs + elif R_SENSE==1: + Mx_cs = CSMETHOD(Mx,r) + My_cs = CSMETHOD(My,r) + Mz_cs = CSMETHOD(Mz,r) + if rank==0: + print('saving the sequences in ' + savepath) + seqname = options['cs']['name'] +'_R' + str(r) + '.npz' + print('sequence name: ' + seqname) + np.savez_compressed( savepath + seqname, x=Mx_cs, y=My_cs,z=Mz_cs) + del Mx_cs,My_cs,Mz_cs + else: + raise Exception('Only implemented for 2-fold SENSE!!') + + + + + + +# THE END + + + diff --git a/codes/Graphics.py b/codes/Graphics.py new file mode 100644 index 0000000..03eb121 --- /dev/null +++ b/codes/Graphics.py @@ -0,0 +1,1581 @@ +import h5py +from dolfin import * +import numpy as np +from matplotlib import pyplot as plt +import scipy as sc +import scipy.io as sio +from mpl_toolkits.mplot3d import Axes3D +from mpl_toolkits.mplot3d import proj3d +import os, sys +import csv +from common import inout +import time + +from matplotlib import rc +#rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']}) +rc('text', usetex=True) + + +def MATmedical(path): + + + datatype = None + + if 'ODV' in path: + datatype = 'odv' + print('The DATA file has the next structure: [x,y,z,coil,CardiacPhase,gradient,NSA] where x \ + and y are the spatial dimensions for the image. z is one because is only one slice. Coil is \ + the number of coils (5), CardiacPhase is the time dimension. Gradient is the number of gradients used in the \ + adquisition (2). NSA is the number of signal average.') + + if 'Phantom' in path: + datatype = 'cib_phantom' + #print('CIB phantom data which contains magnitude and velocity measurements.') + + + if not datatype=='cib_phantom': + STRUC = sio.loadmat(path) + DAT0 = STRUC['data'] + if datatype=='odv': + KDAT0 = STRUC['kdata'] + [row,col,coil,numt2,p4,p5] = DAT0.shape + else: + KDAT0 = STRUC['k_data'] + [row,col,slicesel,coil,numt2,p4,p5] = DAT0.shape + + + z0 = 0 + NSA = 0 + coilnum = 3 + + DAT1 = np.zeros([row,col,numt2],dtype=complex) + DAT2 = np.zeros([row,col,numt2],dtype=complex) + KDAT1 = np.zeros([row,col,numt2],dtype=complex) + KDAT2 = np.zeros([row,col,numt2],dtype=complex) + + if datatype=='odv': + for k in range(numt2): + DAT1[:,:,k] = DAT0[:,:,coilnum,k,0,NSA] + DAT2[:,:,k] = DAT0[:,:,coilnum,k,1,NSA] + KDAT1[:,:,k] = KDAT0[:,:,coilnum,k,0,NSA] + KDAT2[:,:,k] = KDAT0[:,:,coilnum,k,1,NSA] + + return [DAT1,DAT2,KDAT1,KDAT2] + else: + + UK1 = np.zeros([row,col,numt2],dtype=complex) + UK2 = np.zeros([row,col,numt2],dtype=complex) + + for k in range(numt2): + DAT1[:,:,k] = DAT0[:,:,z0,coilnum,k,0,NSA] + DAT2[:,:,k] = DAT0[:,:,z0,coilnum,k,1,NSA] + KDAT1[:,:,k] = KDAT0[:,:,z0,coilnum,k,0,NSA] + KDAT2[:,:,k] = KDAT0[:,:,z0,coilnum,k,1,NSA] + UK1[:,:,k] = np.fft.fftshift(np.fft.fft2(DAT1[:,:,k])) + UK2[:,:,k] = np.fft.fftshift(np.fft.fft2(DAT2[:,:,k])) + + return [DAT1,DAT2,KDAT1,KDAT2] + + elif datatype=='cib_phantom': + + lp = len(path) + for k in range(lp): + if path[lp-1-k]=='/': + name_var = path[lp-k:lp] + break + + name_var = name_var.replace('.mat','') + + if name_var in ['VENC','Segmented_Aorta','voxel_MR']: + STRUC = sio.loadmat(path) + return STRUC[name_var] + else: + #import h5py + arrays = {} + f = h5py.File(path) + for k, v in f.items(): + arrays[k] = np.array(v) + + return arrays[name_var].transpose((1,2,3,0)) + +def Plot_Inlet(): + ''' To Show the flow given at the inlet in the FEM simulation ''' + + Tf = 0.9 + dt = 0.0005 + t = np.linspace(0,Tf,Tf/dt) + U = 60 + Th = 0.37 + DOLFIN_PI = 3.1416 + y = U*np.sin(DOLFIN_PI*t/Th)*(t<=Th) + (Th0: + plt.ylim(ywin) + plt.xlim(xwin) + ax.set_aspect('auto') + + plt.xticks([]) + plt.yticks([]) + #plt.xlabel('$x$',Fontsize=20) + #plt.ylabel('$y$',Fontsize=20) + plt.title('phase ' + str(k)) + plt.show() + + if recorder: + plt.savefig('results/pictures/frame'+str(k),bbox_inches='tight') + + plt.pause(0.101) + plt.clf() + +def PrintSequence2(A,X,Y,Z,x,y,z): + + plt.ion() + fig = plt.figure(figsize=(7, 7), dpi=100) + ax = fig.add_subplot(1,1,1) + [row,col,dep] = A.shape + + + xi = np.linspace(np.min(x),np.max(x),row) + yi = np.linspace(np.min(y),np.max(y),col) + yi2 = np.min(yi) + np.max(yi) - yi + + from matplotlib.colors import LinearSegmentedColormap + basic_cols = ['#ff4e2b', '#000000', '#1893db'] + my_cmap = LinearSegmentedColormap.from_list('mycmap', basic_cols) + ccmap = mpl.cm.Spectral + + xt = x + yt = y + zt = z + yt2 = (np.min(y) + np.max(y)) - y + Dz = (np.max(z) - np.min(z))/dep + + def shiftedColorMap(cmap, start=0, midpoint=0.5, stop=1.0, name='shifted'): + ''' + Function to offset the "center" of a colormap. Useful for + data with a negative min and positive max and you want the + middle of the colormap's dynamic range to be at zero. + + Input + ----- + cmap : The matplotlib colormap to be altered + start : Offset from lowest point in the colormap's range. + Defaults to 0.0 (no lower offset). Should be between + 0.0 and `midpoint`. + midpoint : The new center of the colormap. Defaults to + 0.5 (no shift). Should be between 0.0 and 1.0. In + general, this should be 1 - vmax / (vmax + abs(vmin)) + For example if your data range from -15.0 to +5.0 and + you want the center of the colormap at 0.0, `midpoint` + should be set to 1 - 5/(5 + 15)) or 0.75 + stop : Offset from highest point in the colormap's range. + Defaults to 1.0 (no upper offset). Should be between + `midpoint` and 1.0. + ''' + + cdict = { + 'red': [], + 'green': [], + 'blue': [], + 'alpha': [] + } + + # regular index to compute the colors + reg_index = np.linspace(start, stop, 257) + + # shifted index to match the data + shift_index = np.hstack([ + np.linspace(0.0, midpoint, 128, endpoint=False), + np.linspace(midpoint, 1.0, 129, endpoint=True) + ]) + + for ri, si in zip(reg_index, shift_index): + r, g, b, a = cmap(ri) + + cdict['red'].append((si, r, r)) + cdict['green'].append((si, g, g)) + cdict['blue'].append((si, b, b)) + cdict['alpha'].append((si, a, a)) + + newcmap = mpl.colors.LinearSegmentedColormap(name, cdict) + plt.register_cmap(cmap=newcmap) + + return newcmap + + ccmap = mpl.cm.Spectral + shifted_ccmap = shiftedColorMap(ccmap, midpoint=0.75, name='shifted') + + for k in range(dep): + + zmin = np.min(z) + k*Dz + zmax = zmin + Dz + xs = np.array([]) + ys = np.array([]) + zs = np.array([]) + + for j in range(z.size): + if z[j]>= zmin and z[j] 1: + if os.path.exists(sys.argv[1]): + inputfile = sys.argv[1] + print('Found input file ' + inputfile) + else: + raise Exception('Command line arg given but input file does not exist:' + ' {}'.format(sys.argv[1])) + else: + raise Exception('An input file is required as argument!') + + + start_time = time.time() + + + options = inout.read_parameters(inputfile) + ROUTINE(options) + end_time = time.time() + CLOCK(start_time,end_time) + + + +# END + + + + + + + diff --git a/codes/MATLAB/createU.m b/codes/MATLAB/createU.m new file mode 100644 index 0000000..67eb34e --- /dev/null +++ b/codes/MATLAB/createU.m @@ -0,0 +1,57 @@ +clear all; close all + +folder_name = uigetdir([],'Load Folder...'); + +data = load(strcat(folder_name,'/data.mat')); +SEG = load(strcat(folder_name,'/SEG.mat')); + +data = data.data; +SEG = SEG.SEG; + + +VENC = data.VENC; +VoxelSize = data.voxel_MR; + +vel_AP = data.MR_PCA_AP; +vel_RL = data.MR_PCA_RL; +vel_FH = data.MR_PCA_FH; + +SEG2 = permute(SEG,[2,3,1]); +SEG2 = SEG2(:,:,:); + + +vel_AP_seg = vel_AP.*SEG2(2:end-1,2:end-1,2:end-1); +vel_RL_seg = vel_RL.*SEG2(2:end-1,2:end-1,2:end-1); +vel_FH_seg = vel_FH.*SEG2(2:end-1,2:end-1,2:end-1); + + + + +u_R1 = [] ; +u_R1.x = vel_FH_seg; +u_R1.y = vel_AP_seg; +u_R1.z = vel_RL_seg; +u_R1.VoxelSize = VoxelSize; +save('/home/yeye/Desktop/u_R1.mat','u_R1'); +disp('data saved') +%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% FIGURES +%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +figure +size_vel = size(vel_FH); +for n=1:size_vel(3) + imshow(squeeze(vel_FH_seg(:,:,n,8)),[-100,100],'InitialMagnification',300); + colormap(gca); + pause(0.1) +end +%% +size_seg2 = size(SEG2); +for n=1:size_seg2(3) + imshow(squeeze(SEG2(:,:,n)),'InitialMagnification',300); + colormap(gca); + pause(0.1) +end + + diff --git a/codes/MATLAB/leo/CREATE_MESH.m b/codes/MATLAB/leo/CREATE_MESH.m new file mode 100755 index 0000000..47bab09 --- /dev/null +++ b/codes/MATLAB/leo/CREATE_MESH.m @@ -0,0 +1,14 @@ +% Program to create a structured mesh using the codes of Leo Sok +clear all; close all + +nodes = load('LEO_files/nodes.txt'); +ux = load('LEO_files/ux.txt') ; +uy = load('LEO_files/uy.txt') ; +uz = load('LEO_files/uz.txt') ; +u = sqrt(ux.^2 + uy.^2 + uz.^2); +resol = load('LEO_files/resol.txt') ; +dx = resol(1); dy = resol(2) ; dz = resol(3); + +nodes_masked = maskFEM(nodes,u); +[N,tets,faces] = meshStructTess(nodes_masked,dx,dy,dz,0,0); +writemesh('/home/yeye/Desktop/leomesh',N,tets,faces) \ No newline at end of file diff --git a/codes/MATLAB/leo/maskFEM.m b/codes/MATLAB/leo/maskFEM.m new file mode 100755 index 0000000..b7dca59 --- /dev/null +++ b/codes/MATLAB/leo/maskFEM.m @@ -0,0 +1,19 @@ +function nodes2 = maskFEM(nodes,vel) + +a = []; +b = []; +c = []; +ind = 1; + +for i=1:length(nodes) +if vel(i)>0 + a(ind) = nodes(i,1); + b(ind) = nodes(i,2); + c(ind) = nodes(i,3); + ind = ind +1; +end +end + +nodes2 = [a', b', c']; + + diff --git a/codes/MATLAB/leo/meshStructTess.m b/codes/MATLAB/leo/meshStructTess.m new file mode 100755 index 0000000..4e23f07 --- /dev/null +++ b/codes/MATLAB/leo/meshStructTess.m @@ -0,0 +1,169 @@ +function [nodes, tets, faces, P] = meshStructTess(nodes, dx, dy, dz, check_mesh, plot_mesh) +%% [nodes, tets, faces] = meshStructTess(nodes, dx, dy, dz, check_mesh, plot_mesh) +% Generate a tessalation from a list of structured nodes. +% input: nodes: n times 3 matrix with on the rows the coordinates of +% the n points in the mesh +% dx, dy, dz: the mesh-size in the directions x, y and z +% check_mesh: if true, then it solves a Poisson problem +% plot_mesh: if true, then it plots the mesh +% output: nodes: m times 3 matrix with on the rows the coordinates of +% the m <= n points in the triangulationedi +% tets: l times 4 matrix with on the rows the tetrahedra +% faces: k times 3 matrix with on the rows the triangles of the +% boundary of the mesh +% P: Transformation matrix from input nodes to output nodes. +% Useful also for transforming node-valued functions on +% the input nodes to node-valued functions on the output +% nodes +% +% The triangulation can be plotted using tetramesh(tets,nodes) + + +% compute the minimum and number of points in each direction +if size(nodes,1) < 4 + error('Triangulation needs at least 4 points') +end +mn = min(nodes); +xmin = mn(1); +ymin = mn(2); +zmin = mn(3); + +mn = max(nodes); +xmax = mn(1); +ymax = mn(2); +zmax = mn(3); + +nx = round((xmax-xmin)/dx +1); +ny = round((ymax-ymin)/dy +1); +nz = round((zmax-zmin)/dz +1); + +Nnodes = size(nodes,1); + + +% Define tensor which consist of nodes indices, used for the creation of +% the tetrahedra + +nodes3d = zeros(nx,ny,nz); % preallocate +for i=1:Nnodes + nodes3d(round((nodes(i,1)-xmin)/dx)+1,round((nodes(i,2)-ymin)/dy)+1,round((nodes(i,3)-zmin)/dz)+1)=i; +end + + +disp('Creating Tetrahedra') + +% create tetrahedral mesh in cube, which we will reuse. +ii = 1; +X = zeros(8,3); +for i=0:1 + for j=0:1 + for k=0:1 + X(ii,:) = [i,j,k]; + ii = ii+1; + end + end +end +cubetet = delaunay(X); + +% Run through the mesh +el = 1; +Tetrahedra = zeros(6*(nnz(nodes3d)),4); % preallocate + +for i=1:nx-1 + for j=1:ny-1 + for k=1:nz-1 + % take [i:i+1,j:j+1,k:k+1] as cube + nod = zeros(1,8); % perallocate + + for l = 1:8 + % nod is vector with node indices of cube + nod(l) = nodes3d(i + X(l,1), j + X(l,2), k + X(l,3)); + end + + if nnz(nod) == 8 % then the cube is inside the mesh + tet = nod(cubetet); + else % then there is at least one point of the cube outside the mesh + Xs = X(logical(nod),:); % take only nodes inside the mesh + nodx = nod(logical(nod)); + if nnz(nod) == 4 % 4 nodes, check if points are coplanar + C = cross(Xs(2,:)-Xs(1,:), Xs(3,:)-Xs(1,:)); + cop = logical(dot(C,Xs(4,:)-Xs(1,:))); + % if cop = 0, then points are coplanar end thus no + % tetrahedra exists. + end + if (nnz(nod)>4) || (nnz(nod) == 4 && cop) + % create tetrahedra + tet1 = delaunay(Xs); + tet = nodx(tet1); + else % no tetrahedra exists + tet = []; + end + end + + % add new tetrahedra to list + Tetrahedra(el:el+size(tet,1)-1,:) = tet; + el = el+size(tet,1); + end + end +end + +tets = Tetrahedra(1:el-1,:); % Delete extra preallocated rows. +clear Tetrahedra + +disp([num2str(size(tets,1)), ' tetrahedra created']) + +% Delete nodes which are not in any tetrahedra. +disp('Update mesh') +contr = zeros(size(nodes,1),1); +for i=1:size(tets,1) + for j=1:4 + contr(tets(i,j))=1; + end +end + +nodes = nodes(logical(contr),:); + +% compute P +P = speye(Nnodes); +P = P(logical(contr),:); + +disp([num2str(nnz(~contr)), ' unused nodes in triangulation deleted.']) + +disp('Update tetrahedra') + +% make tetrahedra compatible with new node indices +cumcon = cumsum(~contr)'; +tets = tets - cumcon(tets); + +% create triangles +if size(tets,1) == 0 + warning('No tetrahedra created') + faces = zeros(0,3); +else + disp('Create Triangles') + faces = freeBoundary(triangulation(tets,nodes)); + disp([num2str(size(faces,1)), ' triangles created']) +end + +% checking the mesh by solving a Poisson problem +if check_mesh + % Builds the P1 stiffness matrix from tets and nodes + [A,volumes]=stifness_matrixP1_3D(tets,nodes); + % Check if element volumes may be negative + if any(volumes<=0) + warning('Some elements have zero or negative volume') + end + % solve the Poisson problem with Dirichlet BC + A(2:end,2:end)\ones(size(A(2:end,2:end),1),1); + disp('If there are no warnings, it probably means that the mesh is fine') +end + +% Plots mesh +if plot_mesh + tetramesh(tets,nodes) + xlabel('x') + ylabel('y') + zlabel('z') +end + +end + diff --git a/codes/MATLAB/leo/writemesh.m b/codes/MATLAB/leo/writemesh.m new file mode 100755 index 0000000..a79ad0a --- /dev/null +++ b/codes/MATLAB/leo/writemesh.m @@ -0,0 +1,97 @@ +function writemesh(varargin) +%% writemesh(path, mesh) +% Save triangulation as path.xml and path.msh +% mesh is a struct with fields Pts, Tet, Tri +% alernatively one can use writemesh(path, Pts, Tet, Tri) +% Pts should by a n times 3 matrix consisting points of the mesh +% Tet is the m times 4 matrix consisting the tetrahedra +% Tri is the l times 3 matrix consisting the triangles at the boundary + +if nargin > 3 + mesh.Pts=varargin{2}; + mesh.Tet=varargin{3}; + mesh.Tri=varargin{4}; + writemesh(varargin{1},mesh,varargin(nargin)); + +elseif isstruct(varargin{2}) + rootMeshFile = varargin{1}; + + % NEW FILE + obj = [rootMeshFile,'.msh']; + meshfile = fopen(obj,'w'); + + obj2 = [rootMeshFile,'.xml']; + xmlfile = fopen(obj2,'w'); + + % MESH + fprintf(meshfile,['$MeshFormat','\n']); + fprintf(meshfile,['2.2 0 8','\n']); + fprintf(meshfile,['$EndMeshFormat','\n']); + + fprintf(xmlfile,['','\n']); + fprintf(xmlfile,'\n'); + fprintf(xmlfile,['','\n']); + + mesh = varargin{2}; + + Nodes = mesh.('Pts'); + mesh = rmfield(mesh,'Pts'); + + Nodes = [(1:size(Nodes,1))' Nodes(:,1:3)]; + + % POINTS + if ~strcmp(varargin{nargin},'mute') + disp('Write Points') + end + fprintf(meshfile,['$Nodes','\n']); + fprintf(meshfile,['%i','\n'],size(Nodes,1)); + fprintf(xmlfile,[' ','\n']); + fprintf(xmlfile,[' ','\n'],size(Nodes,1)); + + + fprintf(meshfile,'%i %13.6f %13.6f %13.6f\n',Nodes'); + + Nodes(:,1) = Nodes(:,1) - 1; + + fprintf(xmlfile,' \n',Nodes'); + + fprintf(meshfile,['$EndNodes','\n']); + fprintf(meshfile,['$Elements','\n']); + fprintf(meshfile,['%i','\n'],size(mesh.Tet,1)+size(mesh.Tri,1)); + fprintf(xmlfile,[' ','\n']); + fprintf(xmlfile,[' ','\n'],size(mesh.Tet,1)); + + % Triangles + + if ~strcmp(varargin{nargin},'mute') + disp('Write Triangles') + end + + tri = mesh.('Tri'); + tri = [(1:size(tri,1))' 2*ones(size(tri,1),1) 2*ones(size(tri,1),1) zeros(size(tri,1),1) 2*ones(size(tri,1),1) tri(:,1:3)]; + fprintf(meshfile,'%i %i %i %i %i %i %i %i\n',tri'); + + + + % Tetrahedra + if ~strcmp(varargin{nargin},'mute') + disp('Write Tetrahedra') + end + + tet = mesh.('Tet'); + tet = [(size(tri,1)+1:size(tri,1)+size(tet,1))' 4*ones(size(tet,1),1) 2*ones(size(tet,1),1) zeros(size(tet,1),1) ones(size(tet,1),1) tet(:,1:4)]; + fprintf(meshfile,'%i %i %i %i %i %i %i %i %i\n',tet'); + + tet = mesh.('Tet'); + tet = [(0:size(tet,1)-1)' (tet(:,1:4)-1)]; + fprintf(xmlfile,' \n',tet'); + + + + fprintf(meshfile,['$EndElements','\n']); + fprintf(xmlfile,' \n \n\n'); + + fclose('all'); +end + + diff --git a/codes/MATLAB/load_dicom.m b/codes/MATLAB/load_dicom.m new file mode 100644 index 0000000..35eb488 --- /dev/null +++ b/codes/MATLAB/load_dicom.m @@ -0,0 +1,126 @@ +clear all ; close all +% Load dicom + + +name = 'Ronald' ; + +if strcmp(name, 'Ronald') + path_all = [ + '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Ronald/FH/DICOM/IM_0001', + '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Ronald/AP/DICOM/IM_0001', + '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Ronald/RL/DICOM/IM_0001' + ] ; +end + +if strcmp(name, 'Jeremias') + path_all = [ + '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Jeremias/FH/DICOM/IM_0001', + '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Jeremias/AP/DICOM/IM_0001', + '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190909_Jeremias/RL/DICOM/IM_0001' + ] ; +end + +if strcmp(name, 'Hugo') + path_all = [ + '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Hugo/Dicom/DICOM/IM_0013', + '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Hugo/Dicom/DICOM/IM_0009', + '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Hugo/Dicom/DICOM/IM_0005' + ] ; +end + +for i=1:3 + + if i==1 + %path = '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Paloma/Dicom/DICOM/IM_0013' + disp('Reading the FH component from ...') + path = path_all(1,:) + end + + if i==2 + %path = '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Paloma/Dicom/DICOM/IM_0009' ; + disp('Reading the AP component from ...') + path = path_all(2,:) + end + + if i==3 + %path = '/home/yeye/Desktop/PhD/MEDICAL_DATA/DatosSEPT2019/20190924_Paloma/Dicom/DICOM/IM_0005' ; + disp('Reading the RL component from ...') + path = path_all(3,:) + end + + +I_info = dicominfo(path); +I = double(dicomread(path)); +VENC = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_1.MRVelocityEncodingSequence.Item_1.VelocityEncodingMaximumValue']) ; +heart_rate = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_1.Private_2005_140f.Item_1.HeartRate']); + + +MAG = zeros(size(I,1),size(I,2),I_info.Private_2001_1018,I_info.Private_2001_1017); +PHASE = zeros(size(I,1),size(I,2),I_info.Private_2001_1018,I_info.Private_2001_1017); + +for n=1:size(I,4) + + RI = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.RescaleIntercept']); % intercept + RS = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.RescaleSlope']); % slope + cp = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.Private_2001_1008']); %cp + slc = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.Private_2001_100a']); %scl + id = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_',num2str(n),'.Private_2005_140f.Item_1.Private_2005_106e']); % PCA o FFE + + if strcmp(id,'FFE')==1 + MAG(:,:,slc,cp) = I(:,:,1,n)*RS + RI; + else + PHASE(:,:,slc,cp) = I(:,:,1,n)*RS + RI; + end + +end + + +MASK = double(abs((PHASE==PHASE(1,1,1,1))-1)); +PHASE = PHASE.*MASK; + + +if i==1 + MR_FFE_FH = MAG; + MR_PCA_FH = VENC*PHASE/pi/100; +end + +if i==2 + MR_FFE_AP = MAG; + MR_PCA_AP = VENC*PHASE/pi/100; +end +if i==3 + MR_FFE_RL = MAG; + MR_PCA_RL = VENC*PHASE/pi/100; +end + + +end + + +disp('Saving the data ...') + +spaceslices = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_1.PixelMeasuresSequence.Item_1.SpacingBetweenSlices']); +pixelspacing = eval(['I_info.PerFrameFunctionalGroupsSequence.Item_1.PixelMeasuresSequence.Item_1.PixelSpacing']); + +disp('voxel-size recognized:') +voxel_MR = [pixelspacing(1),pixelspacing(1),spaceslices] + + +data = []; +data.MR_FFE_AP = MR_FFE_AP; +data.MR_FFE_RL = MR_FFE_RL; +data.MR_FFE_FH = MR_FFE_FH; +data.MR_PCA_AP = MR_PCA_AP; +data.MR_PCA_RL = MR_PCA_RL; +data.MR_PCA_FH = MR_PCA_FH; +data.type = 'DAT'; +data.VENC = VENC ; +data.voxel_MR = voxel_MR; +data.heart_rate = heart_rate; + +save('/home/yeye/Desktop/data.mat','data','-v7.3'); +disp('data saved') + + + + \ No newline at end of file diff --git a/codes/MRI.py b/codes/MRI.py new file mode 100644 index 0000000..8ebba11 --- /dev/null +++ b/codes/MRI.py @@ -0,0 +1,1936 @@ +################################################################ +# +# Workspace for MRI analysis of the 4Dflow data +# +# written by Jeremias Garay L: j.e.garay.labra@rug.nl +# Fernanda te amo +# for autoreload in ipython3 +# %load_ext autoreload +# %autoreload 2 +################################################################ +import h5py +from dPdirectestim.dPdirectestim import * +from dolfin import * +import dolfin +import numpy as np +import sys, os +from mpi4py import MPI +from common import inout +import time + +if '2017' in dolfin.__version__: + #class MPI(MPI): + # comm_world = mpi_comm_world() + set_log_active(False) +else: + PRINT_BARS = False + parameters["std_out_all_processes"] = False + +def Etcetera(mode): + if mode=='happy': + print(r""" + /\ /\ + //\\_//\\ ____ + \_ _/ / / + / ^ ^ \ /^^^] + \_\O/_/ [ ] + / \_ [ / + \ \_ / / + [ [ / \/ _/ + _[ [ \ /_/ + """) + + if mode=='serious': + print(r""" + |\__/| + / \ + /_.- -,_\ + \@/ + """) + +def dealiased(VENC,vx,vy,vz): + + vx2 = np.zeros(vx.shape) + vy2 = np.zeros(vy.shape) + vz2 = np.zeros(vz.shape) + + badx_pos = np.where(vx>VENC) + badx_neg = np.where(vx<-VENC) + bady_pos = np.where(vy>VENC) + bady_neg = np.where(vy<-VENC) + badz_pos = np.where(vz>VENC) + badz_neg = np.where(vz<-VENC) + # More than VENC + vx2[badx_pos] = -(2*VENC - vx[badx_pos]) + vy2[bady_pos] = -(2*VENC - vy[bady_pos]) + vz2[badz_pos] = -(2*VENC - vz[badz_pos]) + # Less than VENC + vx2[badx_neg] = 2*VENC + vx[badx_neg] + vy2[bady_neg] = 2*VENC + vy[bady_neg] + vz2[badz_neg] = 2*VENC + vz[badz_neg] + # The rest + otherx = np.where(vx2==0) + othery = np.where(vy2==0) + otherz = np.where(vz2==0) + vx2[otherx] = vx[otherx] + vy2[othery] = vy[othery] + vz2[otherz] = vz[otherz] + + return [vx2,vy2,vz2] + +def checkpoint_norm(MESH,fieldpath,mode): + V = MESH['FEM'] + v = Function(V) + path = fieldpath + 'checkpoint/' + # Checkpoints folders + unsort_indexes = os.listdir(path) + indexes = [int(x) for x in unsort_indexes] + indexes.sort() + norm = np.zeros([len(indexes)]) + j = 0 + comm = MPI.COMM_WORLD + + if mode=='ct': + # CT DATA + for l in indexes: + path1 = path + str(l) + '/u.h5' + hdf = HDF5File(MESH['mesh'].mpi_comm(),path1,'r') + hdf.read(v, 'u/vector_0') + hdf.close() + norm[j] = np.linalg.norm(v.vector().get_local()) + j = j+1 + + if mode=='measurement': + # MEASUREMENT + path2 = fieldpath + 'measurements/' + j = 0 + for l in indexes: + meapath = path2 + 'u' + str(l) + '.h5' + comm = MPI.COMM_WORLD + hdf = HDF5File(MESH['mesh'].mpi_comm(),meapath,'r') + hdf.read(v, 'u/vector_0') + hdf.close() + norm[j] = np.linalg.norm(v.vector().get_local()) + j = j+1 + + + if mode=='roukf': + indexes2 = indexes[1:-1] + j = 1 + for l in indexes2: + path3 = fieldpath + 'roukf/checkpoint/' + str(l) + '/X0.h5' + hdf = HDF5File(MESH['mesh'].mpi_comm(),path3,'r') + hdf.read(v, 'X') + hdf.close() + norm[j] = np.linalg.norm(v.vector().get_local()) + j = j+1 + + + return norm + +def LOADmesh(meshtype,resol=None,boxsize=None,meshpath=None, ranges=None): + + if meshtype=='box_zeros': + # The size of the mesh is given + # The spacements are given + # The x range goes from 0 + [Lx, Ly, Lz] = boxsize + [ddx, ddy, ddz] = resol + x0 = 0 + y0 = 0 + z0 = 0 + x1 = (Lx-1)*ddx + y1 = (Ly-1)*ddy + z1 = (Lz-1)*ddz + mesh0 = BoxMesh(Point(x0, y0, z0), Point(x1, y1, z1), Lx-1, Ly-1, Lz-1) + P1 = FunctionSpace(mesh0, "Lagrange", 1) + # This need to be saved from Leo's construction + Xt = P1.tabulate_dof_coordinates() + #boundaries = MeshFunction('size_t', mesh0 , mesh0.topology().dim() - 1) + #XDMFFile(path2).write(boundaries) + if '2017' in dolfin.__version__: + Xt = np.reshape(Xt, (round(Xt.size/3), 3)) + SqXt = np.zeros(Xt.shape) + for k in range(Xt.shape[0]): + xk = round((Xt[k][0]-x0)/ddx) + yk = round((Xt[k][1]-y0)/ddy) + zk = round((Xt[k][2]-z0)/ddz) + SqXt[k, 0] = int(xk) + SqXt[k, 1] = int(yk) + SqXt[k, 2] = int(zk) + if rank == 0: + print('BOX real created: H = ' + str(mesh0.hmin()) + + ' H_max = ' + str(mesh0.hmax())) + + BOX = {} + BOX['mesh'] = mesh0 + BOX['FEM'] = P1 + BOX['nodes'] = Xt + BOX['seq'] = SqXt.astype(int) + BOX['dimension'] = [Lx,Ly,Lz] + return BOX + + if meshtype=='box_ranges': + # The size of the mesh is given + # The X ranges are given + # The spacements are calculated + [Lx,Ly,Lz] = boxsize + x0 = ranges['x'][0] + x1 = ranges['x'][1] + y0 = ranges['y'][0] + y1 = ranges['y'][1] + z0 = ranges['z'][0] + z1 = ranges['z'][1] + ddx = (x1-x0)/(Lx-1) + ddy = (y1-y0)/(Ly-1) + ddz = (z1-z0)/(Lz-1) + mesh0 = BoxMesh(Point(x0, y0, z0), Point(x1, y1, z1), Lx-1, Ly-1, Lz-1) + P1 = FunctionSpace(mesh0, "Lagrange", 1) + Xt = P1.tabulate_dof_coordinates() + if '2017' in dolfin.__version__: + Xt = np.reshape(Xt, (round(Xt.size/3), 3)) + SqXt = np.zeros(Xt.shape) + for k in range(Xt.shape[0]): + xk = round((Xt[k][0]-x0)/ddx) + yk = round((Xt[k][1]-y0)/ddy) + zk = round((Xt[k][2]-z0)/ddz) + SqXt[k, 0] = int(xk) + SqXt[k, 1] = int(yk) + SqXt[k, 2] = int(zk) + if rank == 0: + print('BOX real created: H = ' + str(mesh0.hmin()) + + ' H_max = ' + str(mesh0.hmax())) + + BOX = {} + BOX['mesh'] = mesh0 + BOX['FEM'] = P1 + BOX['nodes'] = Xt + BOX['seq'] = SqXt.astype(int) + BOX['dimension'] = [Lx, Ly, Lz] + return BOX + + if meshtype=='fix_resolution': + # The size of the mesh is given + # The spacements are given + # The X ranges are given + # The X ranges are redefined from the oldest + [Lx,Ly,Lz] = boxsize + [ddx,ddy,ddz] = resol + x0 = ranges['x'][0] + x1 = ranges['x'][1] + y0 = ranges['y'][0] + y1 = ranges['y'][1] + z0 = ranges['z'][0] + z1 = ranges['z'][1] + deltax = (ddx*(Lx-1) - (x1 - x0))/2 + deltay = (ddy*(Ly-1) - (y1 - y0))/2 + deltaz = (ddz*(Lz-1) - (z1 - z0))/2 + x0_new = x0 - deltax + x1_new = x1 + deltax + y0_new = y0 - deltay + y1_new = y1 + deltay + z0_new = z0 - deltaz + z1_new = z1 + deltaz + + mesh0 = BoxMesh(Point(x0_new,y0_new,z0_new),Point(x1_new,y1_new,z1_new),Lx-1,Ly-1,Lz-1) + P1 = FunctionSpace(mesh0, "Lagrange", 1) + Xt = P1.tabulate_dof_coordinates() # This need to be saved from Leo's construction + #boundaries = MeshFunction('size_t', mesh0 , mesh0.topology().dim() - 1) + #XDMFFile('/home/yeye/Desktop/bound.xdmf').write(boundaries) + if '2017' in dolfin.__version__: + Xt = np.reshape(Xt,( round(Xt.size/3) ,3)) + SqXt = np.zeros(Xt.shape) + for k in range(Xt.shape[0]): + xk = round((Xt[k][0]-x0_new)/ddx) + yk = round((Xt[k][1]-y0_new)/ddy) + zk = round((Xt[k][2]-z0_new)/ddz) + SqXt[k, 0] = int(xk) + SqXt[k, 1] = int(yk) + SqXt[k, 2] = int(zk) + + if rank==0: + print('BOX real created: H = ' + str(mesh0.hmin()) + ' H_max = ' + str(mesh0.hmax()) ) + + BOX = {} + BOX['mesh'] = mesh0 + BOX['FEM'] = P1 + BOX['nodes'] = Xt + BOX['seq'] = SqXt.astype(int) + BOX['dimension'] = [Lx, Ly, Lz] + + return BOX + + if meshtype=='cib': + ddx = 2.00032/1000; ddy = 2.00032/1000; ddz = 2.00526316/1000 + x0 = ddx*17 ; y0 = ddy*17; z0 = ddz + Lx = 75; Ly = 85; Lz = 57 + x1 = x0+ddx*(Lx-1); y1 = y0+ddy*(Ly-1); z1 = z0+ddz*(Lz-1) + + + mesh0 = BoxMesh(Point(x0,y0,z0),Point(x1,y1,z1),Lx-1,Ly-1,Lz-1) + P1 = FunctionSpace(mesh0, "Lagrange", 1) + Xt = P1.tabulate_dof_coordinates() + #XDMFFile('mesh_hex.xdmf').write(mesh0) + # MESH HEXA + #mesh0 = BoxMesh.create([Point(x0,y0,z0), Point(x1,y1,z1)], [Lx, Ly, Lz], CellType.Type.hexahedron) + #P1 = FunctionSpace(mesh0, "DG", 0) + #t = P1.tabulate_dof_coordinates() + + if '2017' in dolfin.__version__: + Xt = np.reshape(Xt,( round(Xt.size/3) ,3)) + SqXt = np.zeros(Xt.shape) + + for k in range(Xt.shape[0]): + xk = round((Xt[k][0]-x0)/ddx) + yk = round((Xt[k][1]-y0)/ddy) + zk = round((Xt[k][2]-z0)/ddz) + SqXt[k,0] = int(xk) + SqXt[k,1] = int(yk) + SqXt[k,2] = int(zk) + if rank==0: + print('BOX cib created: H = ' + str(mesh0.hmin()) + ' H_max = ' + str(mesh0.hmax()) ) + + BOX = {} + BOX['mesh'] = mesh0 + BOX['FEM'] = P1 + BOX['nodes'] = Xt + BOX['seq'] = SqXt.astype(int) + BOX['name'] = 'cib' + BOX['dimension'] = [Lx,Ly,Lz] + return BOX + + if meshtype=='cib_RT': + ddx = 1.8/1000; ddy = 1.8/1000; ddz = 1.8/1000 + x0 = ddx*33 ; y0 = ddy*43; z0 = ddz*13 + Lx = 103; Ly = 63; Lz = 47 + x1 = x0+ddx*(Lx-1); y1 = y0+ddy*(Ly-1); z1 = z0+ddz*(Lz-1) + + mesh0 = BoxMesh(Point(x0,y0,z0),Point(x1,y1,z1),Lx-1,Ly-1,Lz-1) + P1 = FunctionSpace(mesh0, "Lagrange", 1) + Xt = P1.tabulate_dof_coordinates() + #XDMFFile('mesh_hex.xdmf').write(mesh0) + # MESH HEXA + #mesh0 = BoxMesh.create([Point(x0,y0,z0), Point(x1,y1,z1)], [Lx, Ly, Lz], CellType.Type.hexahedron) + #P1 = FunctionSpace(mesh0, "DG", 0) + #t = P1.tabulate_dof_coordinates() + + if '2017' in dolfin.__version__: + Xt = np.reshape(Xt,( round(Xt.size/3) ,3)) + SqXt = np.zeros(Xt.shape) + + for k in range(Xt.shape[0]): + xk = round((Xt[k][0]-x0)/ddx) + yk = round((Xt[k][1]-y0)/ddy) + zk = round((Xt[k][2]-z0)/ddz) + SqXt[k,0] = int(xk) + SqXt[k,1] = int(yk) + SqXt[k,2] = int(zk) + if rank==0: + print('BOX cib created: H = ' + str(mesh0.hmin()) + ' H_max = ' + str(mesh0.hmax()) ) + + BOX = {} + BOX['mesh'] = mesh0 + BOX['FEM'] = P1 + BOX['nodes'] = Xt + BOX['seq'] = SqXt.astype(int) + BOX['name'] = 'cib' + BOX['dimension'] = [Lx,Ly,Lz] + return BOX + + if meshtype=='cib_GM': + ddx = 1.8/1000; ddy = 1.8/1000; ddz = 1.8/1000 + x0 = ddx*20 ; y0 = ddy*43; z0 = -ddz*5 + Lx = 110; Ly = 70; Lz = 43 + x1 = x0+ddx*(Lx-1); y1 = y0+ddy*(Ly-1); z1 = z0+ddz*(Lz-1) + + mesh0 = BoxMesh(Point(x0,y0,z0),Point(x1,y1,z1),Lx-1,Ly-1,Lz-1) + P1 = FunctionSpace(mesh0, "Lagrange", 1) + Xt = P1.tabulate_dof_coordinates() + #XDMFFile('mesh_hex.xdmf').write(mesh0) + # MESH HEXA + #mesh0 = BoxMesh.create([Point(x0,y0,z0), Point(x1,y1,z1)], [Lx, Ly, Lz], CellType.Type.hexahedron) + #P1 = FunctionSpace(mesh0, "DG", 0) + #t = P1.tabulate_dof_coordinates() + + if '2017' in dolfin.__version__: + Xt = np.reshape(Xt,( round(Xt.size/3) ,3)) + SqXt = np.zeros(Xt.shape) + + for k in range(Xt.shape[0]): + xk = round((Xt[k][0]-x0)/ddx) + yk = round((Xt[k][1]-y0)/ddy) + zk = round((Xt[k][2]-z0)/ddz) + SqXt[k,0] = int(xk) + SqXt[k,1] = int(yk) + SqXt[k,2] = int(zk) + if rank==0: + print('BOX cib created: H = ' + str(mesh0.hmin()) + ' H_max = ' + str(mesh0.hmax()) ) + + BOX = {} + BOX['mesh'] = mesh0 + BOX['FEM'] = P1 + BOX['nodes'] = Xt + BOX['seq'] = SqXt.astype(int) + BOX['name'] = 'cib' + BOX['dimension'] = [Lx,Ly,Lz] + return BOX + + if meshtype=='aorta': + + if '.xml' in meshpath: + mesh1 = Mesh(meshpath) + Evel1 = VectorElement('Lagrange', mesh1.ufl_cell(), 1) + V1 = FunctionSpace(mesh1, Evel1) + boundaries1 = MeshFunction( + "size_t", mesh1, mesh1.topology().dim()-1) + if '.h5' in meshpath: + mesh1 = Mesh() + hdf1 = HDF5File(mesh1.mpi_comm(), meshpath, 'r') + hdf1.read(mesh1, '/mesh', False) + boundaries1 = MeshFunction( + 'size_t', mesh1, mesh1.topology().dim() - 1) + hdf1.read(boundaries1, '/boundaries') + Evel1 = VectorElement('Lagrange', mesh1.ufl_cell(), 1) + V1 = FunctionSpace(mesh1, Evel1) + + if rank == 0: + print('AORTA created: h_min = ' + str(mesh1.hmin()) + + ' h_max = ' + str(mesh1.hmax())) + + AORTA = {} + AORTA['mesh'] = mesh1 + AORTA['boundaries'] = boundaries1 + AORTA['FEM'] = V1 + AORTA['name'] = 'aorta' + + return AORTA + + if meshtype == 'leo': + + # LEO + mesh = Mesh(meshpath) + Evel = VectorElement('Lagrange', mesh.ufl_cell(), 1) + #Evel = VectorElement('DG',mesh2.ufl_cell(),0) + V = FunctionSpace(mesh, Evel) + boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1) + marked = True + + if marked: + class Outlet(SubDomain): + def inside(self, x, on_boundary): + return on_boundary and between(x[2], (11.18, 11.4)) + + class Inlet(SubDomain): + def inside(self, x, on_boundary): + return on_boundary and between(x[2], (7.56, 7.8)) and between(x[0], (10, 16)) + + outlet = Outlet() + inlet = Inlet() + boundaries.set_all(0) + outlet.mark(boundaries, 3) + inlet.mark(boundaries, 2) + + if rank == 0: + print('LEO created: H = ' + str(mesh.hmin()) + + ' H_max = ' + str(mesh.hmax())) + + LEO = {} + LEO['mesh'] = mesh + LEO['FEM'] = V + LEO['boundaries'] = boundaries + + return LEO + + if meshtype=='other': + + if '.xml' in meshpath: + mesh1 = Mesh(meshpath) + Evel1 = VectorElement('Lagrange', mesh1.ufl_cell(), 1) + V1 = FunctionSpace(mesh1, Evel1) + boundaries1 = MeshFunction( + "size_t", mesh1, mesh1.topology().dim()-1) + if '.h5' in meshpath: + mesh1 = Mesh() + hdf1 = HDF5File(mesh1.mpi_comm(), meshpath, 'r') + hdf1.read(mesh1, '/mesh', False) + boundaries1 = MeshFunction( + 'size_t', mesh1, mesh1.topology().dim() - 1) + hdf1.read(boundaries1, '/boundaries') + Evel1 = VectorElement('Lagrange', mesh1.ufl_cell(), 1) + V1 = FunctionSpace(mesh1, Evel1) + + if rank == 0: + print('MESH created: h_min = ' + str(mesh1.hmin()) + + ' h_max = ' + str(mesh1.hmax())) + + MESH = {} + MESH['mesh'] = mesh1 + MESH['boundaries'] = boundaries1 + MESH['FEM'] = V1 + + return MESH + +def CREATEsequences(BOX,AORTA,options): + # Checkpoints folders + checkpoint_path = options['sequence']['checkpoint_path']+'checkpoint/' + unsort_indexes = os.listdir(checkpoint_path) + indexes = [int(x) for x in unsort_indexes] + indexes.sort() + [Lx, Ly, Lz] = BOX['dimension'] + + Lt = len(indexes) + Sqx = np.zeros([Lx, Ly, Lz, Lt]) + Sqy = np.zeros([Lx, Ly, Lz, Lt]) + Sqz = np.zeros([Lx, Ly, Lz, Lt]) + + v1x = Function(BOX['FEM']) + v1y = Function(BOX['FEM']) + v1z = Function(BOX['FEM']) + v2 = Function(AORTA['FEM']) + N_sp = int(options['sequence']['sampling_rate']) + + if rank==0: + print('Total number of checkpoints folders: ',len(indexes)) + + if options['sequence']['sampling_rate']>1: + if rank==0: + print('Applying an undersampling of ' + str(N_sp)) + print('list of indexes',indexes) + + for k in range(0, len(indexes), N_sp): + if rank == 0: + print('READING CHECKPOINT FOLDER: ' + str(indexes[k])) + + hdf = HDF5File(AORTA['mesh'].mpi_comm( + ), checkpoint_path + str(indexes[k])+'/u.h5', 'r') + hdf.read(v2, 'u/vector_0') + hdf.close() + + LagrangeInterpolator.interpolate(v1x, v2.sub(0)) + LagrangeInterpolator.interpolate(v1y, v2.sub(1)) + LagrangeInterpolator.interpolate(v1z, v2.sub(2)) + + v1_gather = v1x.vector().gather( + np.array(range(BOX['FEM'].dim()), 'intc')) + v2_gather = v1y.vector().gather( + np.array(range(BOX['FEM'].dim()), 'intc')) + v3_gather = v1z.vector().gather( + np.array(range(BOX['FEM'].dim()), 'intc')) + + SqXt = BOX['seq'] + #SqXt = SqXt.astype(int) + # To Gather the Sorting-Sequence in parallel with Gatherv + loc_X = SqXt[:, 0] + loc_Y = SqXt[:, 1] + loc_Z = SqXt[:, 2] + sendbuf_X = np.array(loc_X) + sendbuf_Y = np.array(loc_Y) + sendbuf_Z = np.array(loc_Z) + sendcounts_X = np.array(comm.gather(len(sendbuf_X), 0)) + sendcounts_Y = np.array(comm.gather(len(sendbuf_Y), 0)) + sendcounts_Z = np.array(comm.gather(len(sendbuf_Z), 0)) + + if rank == 0: + SStot_X = np.zeros([Lx*Ly*Lz], dtype=int) + SStot_Y = np.zeros([Lx*Ly*Lz], dtype=int) + SStot_Z = np.zeros([Lx*Ly*Lz], dtype=int) + else: + SStot_X = None + SStot_Y = None + SStot_Z = None + + comm.Gatherv(sendbuf=sendbuf_X, recvbuf=( + SStot_X, sendcounts_X), root=0) + comm.Gatherv(sendbuf=sendbuf_Y, recvbuf=( + SStot_Y, sendcounts_Y), root=0) + comm.Gatherv(sendbuf=sendbuf_Z, recvbuf=( + SStot_Z, sendcounts_Z), root=0) + + if rank == 0: + for l in range(Lx*Ly*Lz): + Sqx[SStot_X[l], SStot_Y[l], SStot_Z[l], k] = v1_gather[l] + Sqy[SStot_X[l], SStot_Y[l], SStot_Z[l], k] = v2_gather[l] + Sqz[SStot_X[l], SStot_Y[l], SStot_Z[l], k] = v3_gather[l] + + if options['sequence']['sampling_rate']>1: + Sqx = Sqx[:, :, :, 0::N_sp] + Sqy = Sqy[:, :, :, 0::N_sp] + Sqz = Sqz[:, :, :, 0::N_sp] + + if rank == 0: + print('saving the sequences' + options['sequence']['checkpoint_path']) + savepath = options['sequence']['checkpoint_path'] + options['sequence']['name'] + '.npz' + np.savez_compressed(savepath, x=Sqx, y=Sqy, z=Sqz) + + + return [Sqx,Sqy,Sqz] + +def LOADsequences(loadpath): + # for loading existing sequences alocated in loadpath + if rank==0: + print('{reading} ' + loadpath) + + if '.mat' in loadpath: + import scipy.io as sio + S = sio.loadmat(loadpath) + S = S['u_R1'] + Sqx = S['x'][0][0] + Sqy = S['y'][0][0] + Sqz = S['z'][0][0] + else: + S = np.load(loadpath) + Sqx = S['x'] + Sqy = S['y'] + Sqz = S['z'] + + return [Sqx,Sqy,Sqz] + +def SqtoH5(BOX,MESH,Sqx,Sqy,Sqz): + + Xt = BOX['nodes'] + SqXt = BOX['seq'] + P1 = BOX['FEM'] + V = MESH['FEM'] + velocity = {} + Lt = Sqx.shape[3] + VX = V.sub(0).collapse() + VY = V.sub(1).collapse() + VZ = V.sub(2).collapse() + vv_x = Function(VX) + vv_y = Function(VY) + vv_z = Function(VZ) + + if rank == 0: + print('{SQtoH5} total number of timesteps: ' + str(Lt)) + + nit = 0 + nit_tot = (Lt-1)/20 + + for t in range(Lt): + if rank == 0: + if PRINT_BARS: + if np.mod(t, nit_tot) < 1: + sys.stdout.write('\r') + # Progress bar + sys.stdout.write( + "{SQtoH5} [%-40s] %d%%" % ('=='*nit, 100*t/Lt)) + sys.stdout.flush() + nit = nit + 1 + else: + print('iteration number:', t) + + v1_x = Function(P1) + v1_y = Function(P1) + v1_z = Function(P1) + + values_x = np.zeros(v1_x.vector().get_local().shape) + values_y = np.zeros(v1_y.vector().get_local().shape) + values_z = np.zeros(v1_z.vector().get_local().shape) + S0x = Sqx[:, :, :, t] + S0y = Sqy[:, :, :, t] + S0z = Sqz[:, :, :, t] + + for k in range(Xt.shape[0]): + values_x[k] = S0x[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] + values_y[k] = S0y[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] + values_z[k] = S0z[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] + + v1_x.vector()[:] = values_x + v1_y.vector()[:] = values_y + v1_z.vector()[:] = values_z + + LagrangeInterpolator.interpolate(vv_x, v1_x) + LagrangeInterpolator.interpolate(vv_y, v1_y) + LagrangeInterpolator.interpolate(vv_z, v1_z) + + v = Function(V) + assign(v.sub(0), vv_x) + assign(v.sub(1), vv_y) + assign(v.sub(2), vv_z) + + velocity[t] = v + + + return velocity + +def MASKED(S_REF,S_RAW): + + [Sq0x,Sq0y,Sq0z] = S_REF + [Sq1x,Sq1y,Sq1z] = S_RAW + + if Sq0z.shape != Sq1z.shape: + print('Original',Sq0z.shape) + print('To mask',Sq1z.shape) + + raise Exception('Dimension of masking sequences not match!') + + # masked Sq1 according to Sq0. Both sequences must have the same sizes. + Sq2x = np.zeros(Sq0x.shape) + Sq2y = np.zeros(Sq0y.shape) + Sq2z = np.zeros(Sq0z.shape) + + nonzerox = np.where(Sq0x!=0) + nonzeroy = np.where(Sq0y!=0) + nonzeroz = np.where(Sq0z!=0) + + Sq2x[nonzerox] = Sq1x[nonzerox] + Sq2y[nonzeroy] = Sq1y[nonzeroy] + Sq2z[nonzeroz] = Sq1z[nonzeroz] + + return [Sq2x,Sq2y,Sq2z] + +def READcheckpoint(MESH,mode,ckpath,under_rate,name,zeros=True): + + V = MESH['FEM'] + W = MESH['FEM'].sub(0).collapse() + # Checkpoints folders + unsort_indexes = os.listdir(ckpath) + indexes = [int(x) for x in unsort_indexes] + indexes.sort() + + numt = 0 + nit = 0 + nit_tot = (len(indexes)-1)/20 + + from dolfin import HDF5File + if mode == 'u': + velocity = {} + for k in range(0, len(indexes), under_rate): + path = ckpath + str(indexes[k]) + '/'+name+'.h5' + + if zeros: + if indexes[k] < 10 and indexes[k] > 0: + path = ckpath + '0' + str(indexes[k]) + '/'+name+'.h5' + + if rank == 0: + if PRINT_BARS: + if np.mod(numt, nit_tot) < 1: + sys.stdout.write('\r') + sys.stdout.write( + "{vel checkpoint} [%-40s] %d%%" % ('=='*nit, 100*numt/(len(indexes)-1))) + sys.stdout.flush() + nit = nit + 1 + + #if options['P_measurement']=='p2': + #v = Function(V) + #vmix = Function(V) + #v,p = vmix.split() + #from dolfin import HDF5File + #hdf = HDF5File(AORTA['mesh'].mpi_comm(), path, 'r') + #hdf.read(v, 'u') + #hdf.close() + #udat = h5py.File(path, 'r') + #v.vector()[:] = udat['u']['vector_0'] + #else: + + v = Function(V) + comm = MPI.COMM_WORLD + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(v, 'u/vector_0') + hdf.close() + velocity[numt] = v + numt = numt + 1 + + return velocity + + if mode == 'p': + pressure = {} + for k in range(0, len(indexes), under_rate): + path = ckpath + str(indexes[k]) + '/'+name+'.h5' + if zeros: + if indexes[k] < 10 and indexes[k] > 0: + path = ckpath + '0' + str(indexes[k]) + '/'+name+'.h5' + + if rank == 0: + if PRINT_BARS: + if np.mod(numt, nit_tot) < 1: + sys.stdout.write('\r') + sys.stdout.write( + "{press checkpoint} [%-40s] %d%%" % ('=='*nit, 100*numt/(len(indexes)-1))) + sys.stdout.flush() + nit = nit + 1 + p = Function(W) + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(p, 'p/vector_0') + hdf.close() + pressure[numt] = p + numt = numt + 1 + + return pressure + +def CHANGEmesh(MESH,field,mode): + # To change between meshes the same result + field2 = {} + if mode=='velocity': + V = MESH['FEM'] + elif mode=='pressure': + V = MESH['FEM'].sub(0).collapse() + else: + raise Exception('You need to assign a mode for the changing!') + for k in range(len(list(field))): + if rank==0: + print('CHANGING MESH: index',k) + v2 = Function(V) + LagrangeInterpolator.interpolate(v2,field[k]) + field2[k] = v2 + + + return field2 + +def GradP(MESH,pressure,options): + + W = MESH['FEM'].sub(0).collapse() + boundaries = MESH['boundaries'] + barye2mmHg = 1/1333.22387415 + + Tf = Constant(options['Tf']) + dt = Constant(options['Tf']/len(list(pressure))) + NT = int(float(Tf)/float(dt)) + + ind = 0 + t = float(dt) + ds = Measure("ds", subdomain_data=boundaries) + ones = interpolate(Constant(1),W) + DELTA_p = np.zeros([NT]) + + p0 = Function(W) + nit = 0 + nit_tot = (NT-1)/20 + + zero_point = None + if options['reference']['fix_const']: + zero_point = options['reference']['zero_point'] + zero_point = np.array(zero_point) + ndim = W.mesh().topology().dim() + zero_point = W.tabulate_dof_coordinates().reshape((-1, ndim))[0] + + #################### THE TIME LOOP ######################### + while ind < NT: + p0.assign(pressure[ind]) + + if rank == 0: + if PRINT_BARS: + if np.mod(ind, nit_tot) < 1: + sys.stdout.write('\r') + # Progress bar + sys.stdout.write("{ref} [%-40s] %f seg" % + ('=='*nit, round(t, 2))) + sys.stdout.flush() + nit = nit + 1 + + if options['reference']['fix_const']: + p0.vector()[:] -= p0(zero_point) + + i_pin = assemble(p0*ds(2)) + i_pout = assemble(p0*ds(3)) + i_s2 = assemble(ones*ds(2)) + i_s3 = assemble(ones*ds(3)) + + DELTA_p[ind] = barye2mmHg*(i_pin/i_s2-i_pout/i_s3) + ind = ind + 1 + t = t + float(dt) + + if options['reference']['save']: + + if rank==0: + np.savetxt(options['reference']['savepath'],DELTA_p) + print('\n Pressure drop save in ' + options['reference']['savepath']) + +def Fluxes(MESH,velocity,options,bnd): + + mesh = MESH['mesh'] + W = MESH['FEM'] + boundaries = MESH['boundaries'] + n = FacetNormal(mesh) + Tf = Constant(options['Tf']) + rho = Constant(options['density']) + dt = Constant(options['Tf']/len(list(velocity))) + NT = int(float(Tf)/float(dt)) + ind = 0 + t = float(dt) + ds = Measure("ds", subdomain_data=boundaries) + QQ = {} + for k in bnd: + QQ[k] = np.zeros([NT]) + + u = Function(W) + nit = 0 + nit_tot = (NT-1)/20 + #################### THE TIME LOOP ######################### + while indVENC) + INDy = np.any(np.abs(vy)>VENC) + INDz = np.any(np.abs(vz)>VENC) + + if INDx or INDy or INDz: + if not options['phase_contrast']['dealiased']: + if rank==0: + raise Exception('There is some velocity greater or lower than VENC. You should turn dealised true!') + + + if options['phase_contrast']['dealiased']: + [vx_d,vy_d,vz_d] = dealiased(VENC,vx,vy,vz) + else: + vx_d = vx; vy_d = vy; vz_d = vz + + if rank==0: + print('saving velocity as ' + options['phase_contrast']['output']) + np.savez_compressed(options['phase_contrast']['output'], x=vx_d, y=vy_d, z=vz_d) + + if 'resize' in options: + if options['resize']['apply']: + if rank == 0: + print('Changing the size of the sequence') + + [Nx, Ny, Nz] = options['resize']['dim'] + [Sqx, Sqy, Sqz] = LOADsequences(options['resize']['loadseq']) + [Nx0, Ny0, Nz0, Nt] = Sqx.shape + + Sqx_new = np.zeros([Nx, Ny, Nz, Nt]) + Sqy_new = np.zeros([Nx, Ny, Nz, Nt]) + Sqz_new = np.zeros([Nx, Ny, Nz, Nt]) + + ax = int((Nx - Nx0)/2) + bx = ax + Nx0 + ay = int((Ny - Ny0)/2) + by = ay + Ny0 + az = int((Nz - Nz0)/2) + bz = az + Nz0 + + for t in range(Nt): + Sqx_new[ax:bx, ay:by, az:bz, t] = Sqx[:, :, :, t] + Sqy_new[ax:bx, ay:by, az:bz, t] = Sqy[:, :, :, t] + Sqz_new[ax:bx, ay:by, az:bz, t] = Sqz[:, :, :, t] + + if options['resize']['save']: + if rank == 0: + print('saving the sequences' + + options['resize']['savepath']) + np.savez_compressed( + options['resize']['savepath'], x=Sqx_new, y=Sqy_new, z=Sqz_new) + + if 'magnetization_CIB' in options: + if options['magnetization_CIB']['apply']: + if rank == 0: + print('-- Generating k-space for CIB phantom data-- ') + print('Using 0-G gradients') + + import Graphics + masterpath = options['magnetization_CIB']['loadpath'] + mask = Graphics.MATmedical(masterpath+'Segmented_Aorta.mat') + mask = mask.transpose((2, 0, 1)) + VENC = Graphics.MATmedical(masterpath+'VENC.mat') + VENC = VENC[0][0] + if rank == 0: + print('Recognized VENC = ' + str(VENC) + ' cm/s') + + # Magnitudes + magnitude_AP = Graphics.MATmedical(masterpath+'MR_FFE_AP.mat') + magnitude_FH = Graphics.MATmedical(masterpath+'MR_FFE_FH.mat') + magnitude_RL = Graphics.MATmedical(masterpath+'MR_FFE_RL.mat') + # Velocities + vel_AP = Graphics.MATmedical(masterpath+'MR_PCA_AP.mat') + vel_FH = Graphics.MATmedical(masterpath+'MR_PCA_FH.mat') + vel_RL = Graphics.MATmedical(masterpath+'MR_PCA_RL.mat') + + # Cropping images + a_crop = 48 + b_crop = 178 + mask = mask[:, :, a_crop:b_crop] + magnitude_AP = magnitude_AP[:, :, a_crop:b_crop, :] + vel_AP = vel_AP[:, :, a_crop:b_crop, :] + magnitude_RL = magnitude_RL[:, :, a_crop:b_crop, :] + vel_RL = vel_RL[:, :, a_crop:b_crop, :] + magnitude_FH = magnitude_FH[:, :, a_crop:b_crop, :] + vel_FH = vel_FH[:, :, a_crop:b_crop, :] + + # Reference phase + [Nx, Ny, Nz, Nt] = vel_AP.shape + gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei + B0 = 1.5 # Tesla Magnetic Field Strenght + TE = 5e-3 # Echo-time + [X, Y, Z] = np.meshgrid(np.linspace( + 0, Ny, Ny), np.linspace(0, Nx, Nx), np.linspace(0, Nz, Nz)) + frecX = 0.21 + p_ref = np.zeros([Nx, Ny, Nz]) + p_flow_AP = np.zeros([Nx, Ny, Nz]) + p_flow_RL = np.zeros([Nx, Ny, Nz]) + p_flow_FH = np.zeros([Nx, Ny, Nz]) + #Sx2G = np.zeros([Nx,Ny,Nz,Nt],dtype=complex) + #Sy2G = np.zeros([Nx,Ny,Nz,Nt],dtype=complex) + #Sz2G = np.zeros([Nx,Ny,Nz,Nt],dtype=complex) + Sx20 = np.zeros([Nx, Ny, Nz, Nt], dtype=complex) + Sy20 = np.zeros([Nx, Ny, Nz, Nt], dtype=complex) + Sz20 = np.zeros([Nx, Ny, Nz, Nt], dtype=complex) + + treshold = np.mean(magnitude_AP) + #p_ref = (gamma*B0*TE+frecX*X)*(np.abs(magnitude_AP[:,:,:,7])>treshold) + p_ref = (gamma*B0*TE+frecX*X)*(mask > 0.5) + + for t in range(Nt): + p_flow_AP = p_ref + np.pi*vel_AP[:, :, :, t]/VENC + p_flow_RL = p_ref + np.pi*vel_RL[:, :, :, t]/VENC + p_flow_FH = p_ref + np.pi*vel_FH[:, :, :, t]/VENC + + Sx20[:, :, :, t] = magnitude_AP[:, :, :, t] * \ + np.cos(p_ref) + 1j*magnitude_AP[:, :, :, t]*np.sin(p_ref) + Sy20[:, :, :, t] = magnitude_RL[:, :, :, t] * \ + np.cos(p_ref) + 1j*magnitude_RL[:, :, :, t]*np.sin(p_ref) + Sz20[:, :, :, t] = magnitude_FH[:, :, :, t] * \ + np.cos(p_ref) + 1j*magnitude_FH[:, :, :, t]*np.sin(p_ref) + + #Sx2G[:,:,:,t] = magnitude_AP[:,:,:,t]*np.cos(p_flow_AP) + 1j*magnitude_AP[:,:,:,t]*np.sin(p_flow_AP) + #Sy2G[:,:,:,t] = magnitude_RL[:,:,:,t]*np.cos(p_flow_RL) + 1j*magnitude_RL[:,:,:,t]*np.sin(p_flow_RL) + #Sz2G[:,:,:,t] = magnitude_FH[:,:,:,t]*np.cos(p_flow_FH) + 1j*magnitude_FH[:,:,:,t]*np.sin(p_flow_FH) + + if rank == 0: + print('saving the sequences' + + options['magnetization_CIB']['savepath']) + savepath = options['magnetization_CIB']['savepath'] + basepath = savepath.replace('.npz', '') + savepath0 = basepath + '0.npz' + savepathG = basepath + 'G.npz' + np.savez_compressed(savepath0, x=Sx20, y=Sy20, z=Sz20) + #np.savez_compressed(savepathG, x=Sx2G,y=Sy2G,z=Sz2G) + + if 'magnetization' in options: + if options['magnetization']['apply']: + if rank==0: + print('-- Generating k-space -- ') + print('Using 0-G gradients') + + [Sqx,Sqy,Sqz] = LOADsequences(options['magnetization']['loadseq']) + [Nx,Ny,Nz,Nt] = Sqx.shape + VENC = options['magnetization']['VENC'] + gamma = 267.513e6 # rad/Tesla/sec Gyromagnetic ratio for H nuclei + B0 = 1.5 # Tesla Magnetic Field Strenght + TE = 5e-3 # Echo-time + [X,Y,Z] = np.meshgrid(np.linspace(0,Ny,Ny),np.linspace(0,Nx,Nx),np.linspace(0,Nz,Nz)) + frecX = 0.21 # 0.21 without leaps: 0.01 + + Mx = np.zeros(Sqx.shape) + My = np.zeros(Sqx.shape) + Mz = np.zeros(Sqx.shape) + Sx2G = np.zeros(Sqx.shape,dtype=complex) + Sy2G = np.zeros(Sqx.shape,dtype=complex) + Sz2G = np.zeros(Sqx.shape,dtype=complex) + Sx20 = np.zeros(Sqx.shape,dtype=complex) + Sy20 = np.zeros(Sqx.shape,dtype=complex) + Sz20 = np.zeros(Sqx.shape,dtype=complex) + Fx0 = np.zeros(Sqx.shape) + Fy0 = np.zeros(Sqx.shape) + Fz0 = np.zeros(Sqx.shape) + FxG = np.zeros(Sqx.shape) + FyG = np.zeros(Sqx.shape) + FzG = np.zeros(Sqx.shape) + + # Creating M + for t in range(Nt): + Mx[:,:,:,t] = 0.05 + 0.9*(np.abs(Sqx[:,:,:,t])>0.001) + 0.05*np.sqrt(np.abs(Sqx[:,:,:,t])) + My[:,:,:,t] = 0.05 + 0.9*(np.abs(Sqy[:,:,:,t])>0.001) + 0.05*np.sqrt(np.abs(Sqy[:,:,:,t])) + Mz[:,:,:,t] = 0.05 + 0.9*(np.abs(Sqz[:,:,:,t])>0.001) + 0.05*np.sqrt(np.abs(Sqz[:,:,:,t])) + # Organs + Mx[:,:,:,t] = Mx[:,:,:,t] + 1*( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 < 35**2 ) * ( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 > 33**2 ) + Mx[:,:,:,t] = Mx[:,:,:,t] + 0.7*( 0.3*(X-Ny/2+9)**2 + 0.7*(Y-Nx/2-4)**2 < 3**2 )* ( np.abs(Z-Nz/2)<25 ) + My[:,:,:,t] = My[:,:,:,t] + 1*( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 < 35**2 ) * ( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 > 33**2 ) + My[:,:,:,t] = My[:,:,:,t] + 0.7*( 0.3*(X-Ny/2+9)**2 + 0.7*(Y-Nx/2-4)**2 < 3**2 )* ( np.abs(Z-Nz/2)<25 ) + Mz[:,:,:,t] = Mz[:,:,:,t] + 1*( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 < 35**2 ) * ( (X-Ny/2)**2 + 0.75*(Y-Nx/2)**2 > 33**2 ) + Mz[:,:,:,t] = Mz[:,:,:,t] + 0.7*( 0.3*(X-Ny/2+9)**2 + 0.7*(Y-Nx/2-4)**2 < 3**2 )* ( np.abs(Z-Nz/2)<25 ) + FxG[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + np.pi*Sqx[:,:,:,t]/VENC + FyG[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + np.pi*Sqy[:,:,:,t]/VENC + FzG[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + np.pi*Sqz[:,:,:,t]/VENC + Fx0[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + Fy0[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + Fz0[:,:,:,t] = (gamma*B0*TE+frecX*X)*(np.abs(Mz[:,:,:,t])>0.4) + + Sx20[:,:,:,t] = Mx[:,:,:,t]*np.cos(Fx0[:,:,:,t]) + 1j*Mx[:,:,:,t]*np.sin(Fx0[:,:,:,t]) + Sy20[:,:,:,t] = My[:,:,:,t]*np.cos(Fy0[:,:,:,t]) + 1j*My[:,:,:,t]*np.sin(Fy0[:,:,:,t]) + Sz20[:,:,:,t] = Mz[:,:,:,t]*np.cos(Fz0[:,:,:,t]) + 1j*Mz[:,:,:,t]*np.sin(Fz0[:,:,:,t]) + Sx2G[:,:,:,t] = Mx[:,:,:,t]*np.cos(FxG[:,:,:,t]) + 1j*Mx[:,:,:,t]*np.sin(FxG[:,:,:,t]) + Sy2G[:,:,:,t] = My[:,:,:,t]*np.cos(FyG[:,:,:,t]) + 1j*My[:,:,:,t]*np.sin(FyG[:,:,:,t]) + Sz2G[:,:,:,t] = Mz[:,:,:,t]*np.cos(FzG[:,:,:,t]) + 1j*Mz[:,:,:,t]*np.sin(FzG[:,:,:,t]) + + + # Adding noise by cartesian components + for t in range(Nt): + noise1 = np.random.normal(-0.0007,0.049,[Nx,Ny,Nz]) + noise2 = np.random.normal(-0.0007,0.049,[Nx,Ny,Nz]) + Sx20[:,:,:,t] += noise1 + 1j*noise2 + Sy20[:,:,:,t] += noise1 + 1j*noise2 + Sz20[:,:,:,t] += noise1 + 1j*noise2 + + noise1g = np.random.normal(-0.0007,0.049,[Nx,Ny,Nz]) + noise2g = np.random.normal(-0.0007,0.049,[Nx,Ny,Nz]) + Sx2G[:,:,:,t] += noise1g + 1j*noise2g + Sy2G[:,:,:,t] += noise1g + 1j*noise2g + Sz2G[:,:,:,t] += noise1g + 1j*noise2g + + + if rank==0: + print('saving the sequences' + options['magnetization']['savepath']) + savepath = options['magnetization']['savepath'] + basepath = savepath.replace('.npz','') + savepath0 = basepath + '0.npz' + savepathG = basepath + 'G.npz' + np.savez_compressed(savepath0, x=Sx20,y=Sy20,z=Sz20) + np.savez_compressed(savepathG, x=Sx2G,y=Sy2G,z=Sz2G) + + if 'sequence' in options: + if options['sequence']['apply']: + if rank == 0: + print('--- Creating Sequences ---') + + resol = options['sequence']['resol'] + boxsize = options['sequence']['boxsize'] + ranges = options['sequence']['ranges'] + BOX = LOADmesh(options['sequence']['boxtype'], resol=resol, boxsize=boxsize, ranges=ranges) + meshpath = options['sequence']['meshpath'] + AORTA = LOADmesh('aorta', meshpath=meshpath) + [Sqx, Sqy, Sqz] = CREATEsequences(BOX, AORTA, options) + + if 'norms' in options: + if options['norms']['apply']: + [BOX,AORTA,LEO] = CREATEmeshes(options) + if rank==0: + print('Computing Norms between fields') + ct_norm = checkpoint_norm(AORTA,options['norms']['field_path'],'ct') + kal_norm = checkpoint_norm(AORTA,options['norms']['field_path'],'roukf') + + if options['norms']['plot']: + import matplotlib.pyplot as plt + plt.plot(ct_norm, color='blue',marker='.',label='ct') + #plt.plot(meas_norm,label='meas') + plt.plot(kal_norm,color='green',label='kalman') + #plt.ylim([0,8500]) + plt.legend(fontsize=14) + plt.show() + + if 'CIBtoH5' in options: + print('---Converting CIB files into H5---') + if options['CIBtoH5']['apply']: + path_to_cib = options['CIBtoH5']['data_path'] + outpath = options['CIBtoH5']['outpath'] + times = options['CIBtoH5']['times'] + dt = options['CIBtoH5']['dt'] + flip = options['CIBtoH5']['flip'] + + if options['CIBtoH5']['interpolate']: + mesh_path = options['CIBtoH5']['mesh_path'] + CIBtoH5(path_to_cib,times,dt,outpath,interpolate=True,mesh_path=mesh_path,flip=flip) + else: + CIBtoH5(path_to_cib,times,dt,outpath,flip=flip) + + ######################################## + # + # Undersampling + # + ######################################## + if 'cs' in options: + if options['cs']['apply']: + if rank==0: + print('Applying Compressed Sensing') + + + [Sqx, Sqy, Sqz] = LOADsequences(options['cs']['seqpath']) + + + import CS + if options['cs']['short']: + [Mx,My,Mz] = LOADsequences(options['cs']['Mpath']) + CS.undersampling_short(Mx,My,Mz,options) + else: + CS.undersampling(Sqx,Sqy,Sqz,options,options['cs']['savepath']) + + if 'kt-BLAST' in options: + if options['kt-BLAST']['apply']: + if rank==0: + print('Applying kt-BLAST') + import ktBLAST + ktBLAST.undersampling(Sqx,Sqy,Sqz,options,options['kt-BLAST']['savepath']) + + if 'SENSE' in options: + if options['SENSE']['apply']: + if rank==0: + print('-- SENSE reconstruction --') + import SENSE + + [Mx,My,Mz] = LOADsequences(options['SENSE']['Mpath']) + [MxS,MyS,MzS] = SENSE.undersampling(Mx,My,Mz,options) + + if rank==0: + print('saving the sequences' + options['SENSE']['savepath']) + np.savez_compressed(options['SENSE']['savepath'], x=MxS,y=MyS,z=MzS) + + ######################################## + # + # Writing Checkpoint from Sequence + # + ######################################## + + if 'create_checkpoint' in options: + if options['create_checkpoint']['apply']: + print('--- Create Checkpoint ---') + + + mesh_into = options['create_checkpoint']['mesh_into'] + + + if options['create_checkpoint']['extension'] == '.mat': + import scipy.io as sio + seqpath = options['create_checkpoint']['loadseq'] + seqpath = seqpath + options['create_checkpoint']['seqname'] + seqpath = seqpath + '_R1.mat' + S = sio.loadmat(seqpath) + [ddx,ddy,ddz] = S['u_R1']['VoxelSize'][0][0][0] + [ddx,ddy,ddz] = [ddx/10,ddy/10,ddz/10] #mm to cm + S = S['u_R1'] + Sqx = S['x'][0][0] + [Lx,Ly,Lz,Lt] = Sqx.shape + BOX = LOADmesh('box_zeros', resol=[ddx,ddy,ddz], boxsize=[Lx,Ly,Lz]) + del S, Sqx, seqpath + else: + # Box mesh definition + boxsize = options['create_checkpoint']['boxsize'] + resol = options['create_checkpoint']['resol'] + ranges = options['create_checkpoint']['ranges'] + boxtype = options['create_checkpoint']['boxtype'] + if boxtype == 'fix_resolution': + BOX = LOADmesh(boxtype, resol=resol, boxsize=boxsize, ranges=ranges) + + + MESH = LOADmesh('other', meshpath=mesh_into) + + + for r in options['create_checkpoint']['Rseq']: + if rank == 0: + print('Writing checkpoint from sequence ' + + options['create_checkpoint']['seqname'] + '_R' + str(r) + '.npz') + + seqpath = options['create_checkpoint']['loadseq'] + options['create_checkpoint']['seqname'] + \ + '_R' + str(r) + options['create_checkpoint']['extension'] + + [Sx, Sy, Sz] = LOADsequences(seqpath) + + if options['create_checkpoint']['under_rate']>1: + print('Undersampling with factor: ', + options['create_checkpoint']['under_rate']) + Sx = Sx[:, :, :, 0::options['create_checkpoint']['under_rate']] + Sy = Sy[:, :, :, 0::options['create_checkpoint']['under_rate']] + Sz = Sz[:, :, :, 0::options['create_checkpoint']['under_rate']] + + if options['create_checkpoint']['masked']: + [S0x, S0y, S0z] = LOADsequences( + options['create_checkpoint']['masked_seq']) + [Sx, Sy, Sz] = MASKED([S0x, S0y, S0z], [Sx, Sy, Sz]) + del S0x, S0y, S0z + + vel_seq = SqtoH5(BOX, MESH, -Sx, Sy, -Sz) + #vel_seq = SqtoH5(BOX, MESH, Sx, Sy, Sz) + + if rank == 0: + print(' ') + + comm = MESH['mesh'].mpi_comm() + dt = options['create_checkpoint']['dt'] + if options['create_checkpoint']['xdmf']: + xdmf_u = XDMFFile(options['create_checkpoint']['savepath']+'u.xdmf') + + for l in range(len(vel_seq)): + if rank == 0: + print('saving checkpoint', l) + + path = options['create_checkpoint']['savepath'] + \ + 'R' + str(r) + '/checkpoint/{i}/'.format(i=l) + + if options['create_checkpoint']['xdmf']: + vel_seq[l].rename('velocity', 'u') + xdmf_u.write(vel_seq[l], l*dt) + + inout.write_HDF5_data( + comm, path + '/u.h5', vel_seq[l], '/u', t=l*dt) + + ######################################## + # + # Relative Pressure Estimators + # + ######################################## + if 'reference' in options: + if options['reference']['apply']: + if rank == 0: + print('---Applying Reference Pressure Gradient---') + + ckpath = options['reference']['checkpoint_path'] + under_rate = options['reference']['under_rate'] + MESH = LOADmesh('other', meshpath=options['reference']['meshpath']) + zeros = not 'ct' in ckpath + pressure = READcheckpoint( + MESH, 'p', ckpath, under_rate, 'p', zeros=zeros) + GradP(MESH, pressure, options) + + if 'ppe' in options: + + if options['ppe']['apply']: + [BOX,AORTA,LEO] = CREATEmeshes(options) + if rank==0: + print('Applying PPE') + if options['ppe']['mesh_type']=='aorta': + velocity_model = READcheckpoint(AORTA,'v',options,options['checkpoint_path'],'u') + if rank==0: + print(' ') + PPE(AORTA,velocity_model,options,options['ppe']['name']+ '.txt') + if rank==0: + print(' ') + + if options['ppe']['mesh_type']=='leo': + if not options['cs']['apply'] and not options['ppe']['read_check']: + raise Exception('You need first apply a CS sequence') + if options['ppe']['read_check']: + if options['meshname']=='leo': + velocity = READcheckpoint(LEO,'v',options,options['checkpoint_path'],'u') + MESH = LEO + if rank==0: + print(' ') + PPE(MESH,velocity,options,options['ppe']['name']+ '.txt') + if rank==0: + print(' ') + + + if not options['ppe']['read_check']: + for r in options['cs']['R']: + velocity[r] = SqtoH5(BOX,LEO,Sqx_cs[r],Sqy_cs[r],Sqz_cs[r]) + if rank==0: + print(' ') + PPE(LEO,velocity[r],options,options['ppe']['name']+options['ppe']['under_type']+'_R'+str(r)+'.txt') + if rank==0: + print(' ') + + if options['ppe']['mesh_type']=='leoct': + velocity_model = READcheckpoint(AORTA,'v',options,options['checkpoint_path'],'u') + if rank==0: + print(' ') + velocity_model_leo = CHANGEvel(LEO,velocity_model) + if rank==0: + print(' ') + PPE(LEO,velocity_model_leo,options,options['ppe']['name']+ '.txt') + if rank==0: + print(' ') + + if 'ste' in options: + + if options['ste']['apply']: + [BOX,AORTA,LEO] = CREATEmeshes(options) + if rank==0: + print('Applying STE') + + + if options['ste']['mesh_type']=='aorta': + if velocity_model: + STE(AORTA,velocity_model,options,options['ste']['name']+ '.txt') + else: + velocity_model = READcheckpoint(AORTA,'v',options,options['checkpoint_path'],'u') + + if rank==0: + print('') + + + if options['ste']['int']: + STEint(AORTA,velocity_model,options,options['ste']['name']+ '.txt') + else: + STE(AORTA,velocity_model,options,options['ste']['name']+ '.txt') + + if options['ste']['mesh_type']=='leo': + + if not options['cs']['apply'] and not options['ste']['read_check']: + raise Exception('You need first apply a CS sequence') + if options['ste']['read_check']: + if options['meshname']=='leo': + velocity = READcheckpoint(LEO,'v',options,options['checkpoint_path'],'u') + MESH = LEO + if rank==0: + print(' ') + + if options['ste']['int']: + STEint(MESH,velocity,options,options['ste']['name']+ '.txt') + else: + STE(MESH,velocity,options,options['ste']['name']+ '.txt') + if rank==0: + print(' ') + + if not options['ste']['read_check']: + for r in options['cs']['R']: + if velocity[r]: + STE(LEO,velocity[r],options,options['ste']['name']+options['ste']['under_type']+'_R'+str(r)+'.txt') + if rank==0: + print(' ') + else: + velocity[r] = SqtoH5(BOX,LEO,Sqx_cs[r],Sqy_cs[r],Sqz_cs[r]) + if rank==0: + print(' ') + if options['ste']['int']: + STEint(LEO,velocity[r],options,options['ste']['name']+options['ste']['under_type']+'_R'+str(r)+'.txt') + else: + STE(LEO,velocity[r],options,options['ste']['name']+options['ste']['under_type']+'_R'+str(r)+'.txt') + if rank==0: + print(' ') + + if options['ste']['mesh_type']=='leoct': + if velocity_model_leo: + STE(LEO,velocity_model_leo,options,options['ste']['name']+ '.txt') + if rank==0: + print(' ') + else: + velocity_model = READcheckpoint(AORTA,'v',options,options['checkpoint_path'],'u') + if rank==0: + print(' ') + velocity_model_leo = CHANGEvel(LEO,velocity_model) + if rank==0: + print(' ') + if options['ste']['int']: + STEint(LEO,velocity_model_leo,options,options['ste']['name']+ '.txt') + else: + STE(LEO,velocity_model_leo,options,options['ste']['name']+ '.txt') + if rank==0: + print(' ') + + if 'peak_pv' in options: + + if options['peak_pv']['apply']: + import CS + import pickle + import sys + import logging + # DPestim + logging.getLogger().setLevel(logging.INFO) + parameters['form_compiler']['optimize'] = True + parameters['form_compiler']['cpp_optimize'] = True + parameters['form_compiler']['cpp_optimize_flags'] = '-O3 -ffast-math -march=native' + infile_dp = options['peak_pv']['infile_dp'] + estimator = DPDirectEstim(infile_dp) + + barye2mmHg = 1/1333.22387415 + t_star = 0.185 + if rank==0: + print('Computing Velocity and Pressure at Peak Systole') + print('The inlet max occurs at ' + str(t_star) + ' sec') + + + [BOX,AORTA,LEO] = CREATEmeshes(options) + + if options['peak_pv']['mesh_into']=='leo': + MESH = LEO + del AORTA + if options['peak_pv']['mesh_into']=='aorta': + MESH = AORTA + del LEO + + if options['peak_pv']['p_comp']=='error': + if rank==0: + print('Reading Pressure Reference') + P0_PPE = READcheckpoint(MESH,'p',options,options['checkpoint_path'],'p_PPE_impl_stan') + P0_STE = READcheckpoint(MESH,'p',options,options['checkpoint_path'],'p_STE_impl_stan') + + [Sqx,Sqy,Sqz] = LOADsequences(options['peak_pv']['orig_seq']) + Nite = options['peak_pv']['N_CS'] + [row,col,dep,numt] = Sqz.shape + frms_num = 3 + peakslice = options['peak_pv']['peak_slice'] + R = options['peak_pv']['R'] + + v0 = np.sqrt(Sqx[:,:,:,peakslice]**2 + Sqy[:,:,:,peakslice]**2 + Sqz[:,:,:,peakslice]**2) + max0 = np.where(v0==np.max(v0)) + + qqvec = {} + for ss in options['peak_pv']['flux_bnd']: + qqvec[ss] = np.zeros([Nite]) + + ppemax = np.zeros([Nite]) + stemax = np.zeros([Nite]) + vmax = np.zeros([Nite]) + slrdmax = peakslice + # Selecting around the max only + slrdmax = 1 + Sqx = Sqx[:,:,:,peakslice-1:peakslice+2] + Sqy = Sqy[:,:,:,peakslice-1:peakslice+2] + Sqz = Sqz[:,:,:,peakslice-1:peakslice+2] + + for l in range(len(R)): + if rank==0: + print('Peak Systole velocity and pressure at R = ' + str(R[l])) + + sx_cs = np.zeros([row,col,dep,frms_num,Nite]) + sy_cs = np.zeros([row,col,dep,frms_num,Nite]) + sz_cs = np.zeros([row,col,dep,frms_num,Nite]) + + for k in range(Nite): + if rank==0: + print('CS iteration number ' + str(k+1)) + [xcs,ycs,zcs] = CS.undersampling_peakpv(Sqx,Sqy,Sqz,options,R[l]) + sx_cs[:,:,:,:,k] = xcs + sy_cs[:,:,:,:,k] = ycs + sz_cs[:,:,:,:,k] = zcs + vk = np.sqrt(sx_cs[:,:,:,1,k]**2 + sy_cs[:,:,:,1,k]**2 + sz_cs[:,:,:,1,k]**2) + vmax[k] = vk[max0] + + if rank==0: + print('\n CS done') + + # To write the checkpoints + vel_seq = SqtoH5(BOX,MESH,sx_cs[:,:,:,:,k],sy_cs[:,:,:,:,k],sz_cs[:,:,:,:,k]) + comm = MESH['mesh'].mpi_comm() + + # Computing the Fluxes + if rank==0: + print('\n Computing the Flux') + QQ = Fluxes(MESH,vel_seq,options,options['peak_pv']['flux_bnd']) + + for ss in options['peak_pv']['flux_bnd']: + qqvec[ss][k] = QQ[ss][slrdmax] + + if rank==0: + print('\n Writing checkpoints') + + for ns in range(len(vel_seq)): + pathss = options['peak_pv']['savepath'] + 'H5/checkpoint/{i}/'.format(i=ns) + if l<10 and l>0: + pathss = options['peak_pv']['savepath'] + 'H5/checkpoint/0{i}/'.format(i=ns) + inout.write_HDF5_data(comm, pathss + '/u.h5', vel_seq[ns], '/u', t=0) + if rank==0: + print('\n The checkpoints were wrote') + + # Computing the Pressure Drop + estimator.estimate() + # Reading the results + if options['peak_pv']['p_comp']=='peak': + ppe_raw = open(options['peak_pv']['savepath'] + '/H5/pdrop_PPE_impl_stan.dat','rb') + ste_raw = open(options['peak_pv']['savepath'] + '/H5/pdrop_STE_impl_stan.dat','rb') + ppe = pickle.load(ppe_raw)['pdrop']*(-barye2mmHg) + ste = pickle.load(ste_raw)['pdrop']*(-barye2mmHg) + p1max[k] = ppe[slrdmax] + p2max[k] = ste[slrdmax] + elif options['peak_pv']['p_comp']=='error': + PPE = READcheckpoint(MESH,'p',options,options['peak_pv']['savepath']+'H5/checkpoint/','p_PPE_impl_stan') + STE = READcheckpoint(MESH,'p',options,options['peak_pv']['savepath']+'H5/checkpoint/','p_STE_impl_stan') + ppe_vec_0 = P0_PPE[peakslice].vector().get_local() - P0_PPE[peakslice].vector().get_local()[0] + ppe_vec = PPE[slrdmax].vector().get_local() - PPE[slrdmax].vector().get_local()[0] + ste_vec_0 = P0_STE[peakslice].vector().get_local() - P0_STE[peakslice].vector().get_local()[0] + ste_vec = STE[slrdmax].vector().get_local() - STE[slrdmax].vector().get_local()[0] + ppemax[k] = np.linalg.norm(ppe_vec_0 - ppe_vec)/np.linalg.norm(ppe_vec_0) + stemax[k] = np.linalg.norm(ste_vec_0 - ste_vec)/np.linalg.norm(ste_vec_0) + else: + raise Exception('Pressure computation not recognize!') + + + # VELOCITIES + vmean = np.mean(vmax) + vstd = np.std(vmax) + # PRESSURES + ppemean = np.mean(ppemax) + stemean = np.mean(stemax) + ppestd = np.std(ppemax) + stestd = np.std(stemax) + + + if options['peak_pv']['save']: + if rank==0: + print('\n saving the files in ' + options['peak_pv']['savepath']) + for ss in options['peak_pv']['flux_bnd']: + np.savetxt( options['peak_pv']['savepath'] + 'Q'+str(ss) +'_R'+str(R[l])+'.txt', [np.mean(qqvec[ss])]) + np.savetxt( options['peak_pv']['savepath'] + 'Qstd'+str(ss) +'_R'+str(R[l])+'.txt', [np.std(qqvec[ss])]) + np.savetxt( options['peak_pv']['savepath'] + 'ppemean_R'+str(R[l])+'.txt', [ppemean] ) + np.savetxt( options['peak_pv']['savepath'] + 'stemean_R'+str(R[l])+'.txt', [stemean] ) + np.savetxt( options['peak_pv']['savepath'] + 'vmean_R'+str(R[l])+'.txt', [vmean] ) + np.savetxt( options['peak_pv']['savepath'] + 'ppestd_R'+str(R[l])+'.txt', [ppestd] ) + np.savetxt( options['peak_pv']['savepath'] + 'stestd_R'+str(R[l])+'.txt', [stestd] ) + np.savetxt( options['peak_pv']['savepath'] + 'vstd_R'+str(R[l])+'.txt', [vstd] ) + + if 'change_mesh' in options: + if options['change_mesh']['apply']: + if rank == 0: + print('Changing '+options['change_mesh']['mode'] + ' from: ') + print('Mesh input: ' + options['change_mesh']['mesh_in']) + print('Mesh output: ' + options['change_mesh']['mesh_out']) + + MESH_in = LOADmesh( + 'other', meshpath=options['change_mesh']['mesh_in']) + MESH_out = LOADmesh( + 'other', meshpath=options['change_mesh']['mesh_out']) + + ckpath = options['change_mesh']['checkpoint_path'] + under_rate = options['change_mesh']['under_rate'] + mode = options['change_mesh']['mode'] + + origin = READcheckpoint(MESH_in, mode, ckpath, under_rate, mode, zeros=False) + + # To change between meshes the same result + changed = {} + #W = LEO['FEM'].sub(0).collapse() + comm = MESH_out['mesh'].mpi_comm() + v2 = Function(MESH_out['FEM']) + + for k in range(len(list(origin))): + if rank == 0: + print('CHANGING: index', k) + LagrangeInterpolator.interpolate(v2, origin[k]) + changed[k] = v2 + print(max(origin[k].vector().get_local()),max(v2.vector().get_local())) + + + dt = options['change_mesh']['dt'] + if options['create_checkpoint']['xdmf']: + xdmf_u = XDMFFile(options['change_mesh']['savepath']+'u.xdmf') + for l in range(len(changed)): + if rank == 0: + print('saving checkpoint', l) + path = options['change_mesh']['savepath'] + \ + 'R1/checkpoint/{i}/'.format(i=l) + writepath = path + '/'+options['change_mesh']['mode']+'.h5' + inout.write_HDF5_data( + comm, writepath, changed[l] ,'/'+options['change_mesh']['mode'], t=l*dt) + + if options['create_checkpoint']['xdmf']: + changed[l].rename('velocity', 'u') + xdmf_u.write(changed[l], l*dt) + + + +if __name__ == '__main__': + + comm = MPI.COMM_WORLD + size = comm.Get_size() + rank = comm.Get_rank() + + if len(sys.argv) > 1: + if os.path.exists(sys.argv[1]): + inputfile = sys.argv[1] + if rank==0: + print('Found input file ' + inputfile) + else: + Etcetera('serious') + raise Exception('Command line arg given but input file does not exist:' + ' {}'.format(sys.argv[1])) + else: + Etcetera('serious') + raise Exception('An input file is required as argument!') + + + user = 'p283370' # Default's user + if 'Zion' in os.popen('hostname').read(): + user = 'yeye' + np.set_printoptions(threshold=5) + + if 'fwn-bborg-5-166' in os.popen('hostname').read(): + user = 'p283370' + + if rank==0: + print('Welcome user {uss}'.format(uss=user)) + Etcetera('happy') + + + start_time = time.time() + + options = inout.read_parameters(inputfile) + SCANNER(options) + + end_time = time.time() + CLOCK(rank,start_time,end_time) + diff --git a/codes/PostCheck.py b/codes/PostCheck.py new file mode 100644 index 0000000..1674730 --- /dev/null +++ b/codes/PostCheck.py @@ -0,0 +1,1072 @@ +from dolfin import * +import numpy as np +import sys +import os +from mpi4py import MPI +from common import inout + + +def LOADmesh(pathtomesh): + + mesh = Mesh() + hdf = HDF5File(mesh.mpi_comm(), pathtomesh, 'r') + hdf.read(mesh, '/mesh', False) + boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1) + hdf.read(boundaries, '/boundaries') + + #Evel = VectorElement('P',mesh.ufl_cell(),1) + #V = FunctionSpace(mesh,Evel) + + #V1 = VectorElement('P', mesh.ufl_cell(), 1) + #V = FunctionSpace(mesh, V1) + P1 = FiniteElement('P', mesh.ufl_cell(), 1) + V1 = VectorElement('P', mesh.ufl_cell(), 1) + V = FunctionSpace(mesh, V1*P1) + + if rank == 0: + print('MESH ' + ' created: H = ' + + str(mesh.hmin()) + ' H_max = ' + str(mesh.hmax())) + + MESH = {} + MESH['mesh'] = mesh + MESH['FEM'] = V + MESH['boundaries'] = boundaries + + return MESH + +def LOADsequences(loadpath): + # for loading existing sequences alocated in loadpath + if rank == 0: + print('{reading} ' + loadpath) + S = np.load(loadpath) + Sqx = S['x'] + Sqy = S['y'] + Sqz = S['z'] + return [Sqx, Sqy, Sqz] + +def WORKcheck(MESH, mode, output_path, checkpoint_path, filename, outname, options): + + #from dolfin import HDF5File + V = MESH['FEM'] + W = MESH['FEM'].sub(0).collapse() + # Checkpoints folders + unsort_indexes = os.listdir(checkpoint_path) + indexes = [int(x) for x in unsort_indexes] + indexes.sort() + + + if mode == 'u': + v = Function(W) + xdmf_u = XDMFFile(output_path+outname+'.xdmf') + dt = options['Velocity']['dt'] + for k in range(0,len(indexes),options['Velocity']['undersampling']): + te = k*dt + path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' + v.rename('velocity', outname) + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + if 'w' in filename: + hdf.read(v, 'w/vector_0') + else: + hdf.read(v, 'u/vector_0') + hdf.close() + xdmf_u.write(v, te) + + if mode == 'w': + xdmf_u = XDMFFile(output_path+outname+'.xdmf') + for k in range(0, len(indexes), 1): + path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' + v = Function(V) + v.rename('waux', outname) + comm = MPI.COMM_WORLD + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(v, 'w/vector_0') + hdf.close() + wvec = v.vector().get_local() + wnorm = wvec - np.mean(wvec) + v.vector()[:] = wnorm + xdmf_u.write(v, te) + te = te + dt + numt = numt + 1 + + if mode == 'gradw': + xdmf_u = XDMFFile(output_path+outname+'.xdmf') + for k in range(0, len(indexes), 1): + path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' + v = Function(V) + gw = Function(W) + gw.rename(outname, outname) + comm = MPI.COMM_WORLD + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(v, 'w/vector_0') + hdf.close() + wvec = v.vector().get_local() + wnorm = wvec - np.mean(wvec) + v.vector()[:] = wnorm + gw.assign(project(sqrt(inner(grad(v), grad(v))), W)) + xdmf_u.write(gw, te) + te = te + dt + numt = numt + 1 + + if mode == 'p' or mode == 'p_cib': + xdmf_p = XDMFFile(output_path+outname+'.xdmf') + for k in range(0, len(indexes), 1): + path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' + + if filename == 'p': + if k < 10 and k > 0: + path = checkpoint_path + '0' + \ + str(indexes[k]) + '/'+filename+'.h5' + + p = Function(W) + if mode == 'p': + barye2mmHg = 1/1333.22387415 + if mode == 'p_cib': + barye2mmHg = -0.00750062 + p.rename('pressure', outname) + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(p, 'p/vector_0') + hdf.close() + + p1vec = p.vector().get_local() + p1vec = p1vec - np.mean(p1vec) + p.vector()[:] = p1vec*barye2mmHg + xdmf_p.write(p, te) + te = te + dt + numt = numt + 1 + + if mode == 'divu': + xdmf_u = XDMFFile(output_path+outname+'.xdmf') + for k in range(0, len(indexes), 1): + path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' + + if indexes[k] < 10 and indexes[k] > 0: + path = checkpoint_path + '0' + \ + str(indexes[k]) + '/'+filename+'.h5' + + v = Function(V) + dv = Function(W) + + dv.rename('div', outname) + comm = MPI.COMM_WORLD + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(v, 'u/vector_0') + + hdf.close() + dv.assign(project(div(v), W)) + xdmf_u.write(dv, te) + te = te + dt + numt = numt + 1 + +def READcheckpoint(MESH, mode, output_path, checkpoint_path, filename, outname, options, flow=False, bnds=None): + + #from dolfin import HDF5File + V = MESH['FEM'] + W = MESH['FEM'].sub(0).collapse() + # Checkpoints folders + unsort_indexes = os.listdir(checkpoint_path) + indexes = [int(x) for x in unsort_indexes] + indexes.sort() + + if flow: + QQ = {} + ds = Measure('ds', domain=MESH['mesh'], + subdomain_data=MESH['boundaries']) + n = FacetNormal(MESH['mesh']) + for bb in bnds: + QQ[bb] = [] + + + if mode == 'interpolation': + dt_new = options['Temporal-Interpolation']['dt_new'] + dt = options['Temporal-Interpolation']['dt'] + Nfolder = len(indexes) + Tf = dt*(Nfolder-1) + Nfolder_new = np.int(Tf/dt_new + 1) + indexes_new = [k for k in range(1,Nfolder_new+1)] + + u = Function(W) + unew = Function(W) + u.rename('velocity', outname) + unew.rename('velocity', outname) + + if options['Temporal-Interpolation']['xdmf']: + xdmf_u = XDMFFile(output_path+'u.xdmf') + + uvec = {} + # Reading first all the points for every original time-steps + for k in indexes: + path = checkpoint_path + str(k) + '/'+filename+'.h5' + print('Reading timestep number: ' + str(k)) + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + if 'w' in filename: + hdf.read(u, 'w/vector_0') + else: + hdf.read(u, 'u/vector_0') + hdf.close() + uvec[k] = u.vector().get_local() + + times_old = np.linspace(0,Tf,Nfolder) + times_new = np.linspace(0,Tf,Nfolder_new) + velnodes = np.zeros(times_old.size) + velnodes_new = np.zeros(times_new.size) + uvec_new = {} + + for k in indexes_new: + uvec_new[k] = np.zeros(uvec[1].size) + + from scipy.interpolate import interp1d + print('Interpolating in every node across time') + # FOR every single node!!! + for l in range(len(uvec[1])): + for k in range(len(velnodes)): + velnodes[k] = uvec[indexes[k]][l] + inter_f = interp1d(times_old, velnodes , kind='cubic') + velnodes_new = inter_f(times_new) + for k in range(len(velnodes_new)): + uvec_new[indexes_new[k]][l] = velnodes_new[k] + print('Interpolation done') + + for k in indexes_new: + print('Writing timestep number: ' + str(k) ) + + unew.vector()[:] = uvec_new[k][:] + + write_path = output_path + 'checkpoint/{i}/'.format(i=k) + hdf2 = HDF5File(MESH['mesh'].mpi_comm(), write_path + 'u.h5', 'w') + hdf2.write(unew, '/u', float((k-1)*dt_new)) + hdf2.close() + + if options['Temporal-Interpolation']['xdmf']: + xdmf_u.write(unew, (k-1)*dt_new) + + + if mode == 'average': + N_av = options['Temporal-Average']['subsampling_rate'] + dt = options['Temporal-Average']['dt'] + ks = 1 + mean_fill = False + u = Function(W) + usub = Function(W) + u.rename('velocity', outname) + usub.rename('velocity', outname) + if options['Temporal-Average']['xdmf']: + xdmf_u = XDMFFile(output_path+'test.xdmf') + + for k in range(len(indexes)): + path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' + print('Reading timestep number: ' + str(indexes[k])) + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + if 'w' in filename: + hdf.read(u, 'w/vector_0') + else: + hdf.read(u, 'u/vector_0') + hdf.close() + + if not mean_fill: + mean_u = u.vector().get_local()/N_av + mean_fill = True + else: + mean_u += u.vector().get_local()/N_av + + if np.mod(k+1,N_av) == 0: + mean_fill = False + usub.vector()[:] = mean_u + print('saving timestep number: ' , ks) + write_path = output_path + 'checkpoint/{i}/'.format(i=ks) + hdf2 = HDF5File(MESH['mesh'].mpi_comm(), write_path + 'u.h5', 'w') + hdf2.write(usub, '/u', float((k-N_av+1)*dt)) + hdf2.close() + if options['Temporal-Average']['xdmf']: + xdmf_u.write(usub, (k-N_av+1)*dt) + ks +=1 + + if mode == 'u': + xdmf_u = XDMFFile(output_path+outname+'.xdmf') + for k in range(0, len(indexes), 1): + path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' + u = Function(V) + u.rename('velocity', outname) + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(u, 'u/vector_0') + hdf.close() + if flow: + for bb in bnds: + QQ[bb].append(assemble(dot(u, n)*ds(bb))) + + xdmf_u.write(u, k) + te = te + dt + + if flow: + for bb in bnds: + np.savetxt(output_path+'flowrate'+str(bb)+'.txt', QQ[bb]) + + if mode == 'w': + xdmf_u = XDMFFile(output_path+outname+'.xdmf') + for k in range(0, len(indexes), 1): + path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' + v = Function(V) + v.rename('waux', outname) + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(v, 'w/vector_0') + hdf.close() + wvec = v.vector().get_local() + wnorm = wvec - np.mean(wvec) + v.vector()[:] = wnorm + xdmf_u.write(v, k) + + if mode == 'gradw': + xdmf_u = XDMFFile(output_path+outname+'.xdmf') + for k in range(0, len(indexes), 1): + path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' + v = Function(V) + gw = Function(W) + gw.rename(outname, outname) + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(v, 'w/vector_0') + hdf.close() + wvec = v.vector().get_local() + wnorm = wvec - np.mean(wvec) + v.vector()[:] = wnorm + gw.assign(project(sqrt(inner(grad(v), grad(v))), W)) + xdmf_u.write(gw, k) + + if mode == 'p' or mode == 'p_cib': + xdmf_p = XDMFFile(output_path+outname+'.xdmf') + for k in range(0, len(indexes), 1): + path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' + p = Function(W) + if mode == 'p': + barye2mmHg = -1/1333.22387415 + if mode == 'p_cib': + barye2mmHg = -0.00750062 + p.rename('pressure', outname) + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(p, 'p/vector_0') + hdf.close() + p1vec = p.vector().get_local()*barye2mmHg + p.vector()[:] = p1vec + xdmf_p.write(p, k) + + if mode == 'divu': + xdmf_u = XDMFFile(output_path+outname+'.xdmf') + for k in range(0, len(indexes), 1): + path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' + v = Function(V) + dv = Function(W) + dv.rename('div', outname) + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(v, 'u/vector_0') + hdf.close() + dv.assign(project(div(v), W)) + xdmf_u.write(dv, k) + +def ERRORmap(MESH, mode, outpath, reference_path, checkpoint_path, refname,comname,options): + + from dolfin import HDF5File + V = MESH['FEM'] + W = MESH['FEM'].sub(0).collapse() + unsort_indexes = os.listdir(checkpoint_path) + indexes = [int(x) for x in unsort_indexes] + indexes.sort() + unsort_indexes0 = os.listdir(reference_path) + indexes0 = [int(x) for x in unsort_indexes0] + indexes0.sort() + + + if len(indexes)!=len(indexes0): + raise Exception('The lengh of the checkpoints are not the same!') + + if mode =='curves': + + if options['Error-curves']['undersampling']>1: + print('undersampling in the reference measurements!') + n_under = options['Error-curves']['undersampling'] + indexes0 = indexes0[0::n_under] + + u = Function(W) + w = Function(W) + for typs in options['Error-curves']['type']: + ucomp = [] + wcomp = [] + times = [] + dt = options['Error-curves']['dt'] + for k in range(1,len(indexes)): + path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' + path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' + + u.rename('meas', 'meas') + w.rename('w', 'w') + hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') + hdf_w.read(w, 'w/vector_0') + hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') + hdf_u.read(u, 'u/vector_0') + hdf_u.close() + hdf_w.close() + u_vec = u.vector().get_local() + w_vec = w.vector().get_local() + print('computing error in timestep numer',k) + + if typs == 'mean': + ucomp.append(np.mean(abs(u_vec))) + wcomp.append(np.mean(abs(w_vec))) + elif typs == 'max': + ucomp.append(np.max(abs(u_vec))) + wcomp.append(np.max(abs(w_vec))) + else: + raise Exception('No defined type for curve printing!') + + times.append(k*dt) + + np.savetxt(outpath+'u' +typs+'.txt',ucomp) + np.savetxt(outpath+'w' +typs+'.txt',wcomp) + np.savetxt(outpath+'times.txt',times) + + if mode == 'histogram': + from pathlib import Path + import pickle + u = Function(V) + w = Function(V) + errors = {} + + for k in range(len(indexes)): + path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' + path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' + + if indexes0[k] < 10: + path_u = reference_path + '0' + \ + str(indexes0[k]) + '/'+refname+'.h5' + + u.rename('meas', 'meas') + w.rename('w', 'w') + hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') + hdf_w.read(w, 'w/vector_0') + hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') + hdf_u.read(u, 'u/vector_0') + hdf_u.close() + hdf_w.close() + + u_vec = u.vector().get_local() + w_vec = w.vector().get_local() + errors[k] = np.zeros(u_vec.size) + + for l in range(len(errors[k])): + #errors[k][l] = np.nan + if u_vec[l]<1e-9: + errors[k][l] = -1 + else: + eta = np.abs(w_vec[l]/u_vec[l]) + if np.abs(eta)>50: + errors[k][l] = -1 + else: + errors[k][l] = eta + + + write_path = Path(outpath) + fpath = write_path.joinpath('errors.dat') + pickle.dump(errors, fpath.open('wb')) + + if mode == 'h5': + xdmf_u = XDMFFile(outpath+'meas.xdmf') + #xdmf_ucor = XDMFFile(output_path+'ucor.xdmf') + xdmf_w = XDMFFile(outpath+'w.xdmf') + ds = Measure("ds", subdomain_data=MESH['boundaries']) + + u = Function(W) + w = Function(W) + #ucor = Function(W) + + for k in range(0, len(indexes), 1): + path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' + path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' + + if indexes0[k] < 10: + path_u = reference_path + '0' + \ + str(indexes0[k]) + '/'+refname+'.h5' + + u.rename('meas', 'meas') + w.rename('w', 'w') + #ucor.rename('ucor', 'ucor') + hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') + hdf_w.read(w, 'w/vector_0') + hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') + hdf_u.read(u, 'u/vector_0') + + hdf_u.close() + hdf_w.close() + + #u_vec = u.vector().get_local() + #w_vec = w.vector().get_local() + #ucor.vector()[:] = u_vec + w_vec + + xdmf_u.write(u, k) + #xdmf_ucor.write(ucor, k) + xdmf_w.write(w, k) + + if mode == 'colormap': + colormap = XDMFFile(outpath+'colormap.xdmf') + #ds = Measure("ds", subdomain_data=MESH['boundaries']) + u = Function(W) + w = Function(W) + cm = Function(W) + dt = options['Corrector']['dt'] + + for k in range(len(indexes)): + print('making the colormap in the time',np.round(k*dt,2)) + path_w = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' + path_u = reference_path + str(indexes0[k]) + '/'+refname+'.h5' + u.rename('meas', 'meas') + w.rename('w', 'w') + cm.rename('color','color') + hdf_w = HDF5File(MESH['mesh'].mpi_comm(),path_w,'r') + hdf_w.read(w, 'w/vector_0') + hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') + hdf_u.read(u, 'u/vector_0') + hdf_u.close() + hdf_w.close() + uvec = u.vector().get_local() + wvec = w.vector().get_local() + cm.vector()[:] = np.sqrt((uvec - wvec)**2) + colormap.write(cm, k*dt) + + if mode == 'error_u': + xdmf_u = XDMFFile(output_path+'error_u.xdmf') + for k in range(0, len(indexes), 1): + path = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' + path0 = reference_path + str(indexes0[k]) + '/'+refname+'.h5' + + if k < 10: + path = checkpoint_path + '0'+str(indexes[k]) + '/u.h5' + path0 = reference_path + '0'+str(indexes[k]) + '/u.h5' + + u = Function(V) + u0 = Function(V) + + eu = Function(W) + eu.rename('error_u', 'error_u') + + comm = MPI.COMM_WORLD + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(u, 'u/vector_0') + hdf0 = HDF5File(MESH['mesh'].mpi_comm(), path0, 'r') + hdf0.read(u0, 'u/vector_0') + hdf0.close() + hdf.close() + + eu.assign(project(sqrt(inner(u-u0, u-u0)), W)) + + xdmf_u.write(eu, te) + te = te + dt + numt = numt + 1 + + if mode == 'error_p' or mode == 'error_p_cib': + xdmf_p = XDMFFile(output_path+outname+'.xdmf') + for k in range(0, len(indexes0), 1): + path = checkpoint_path + str(indexes[k]) + '/'+comname+'.h5' + path0 = reference_path + str(indexes0[k]) + '/'+refname+'.h5' + + if refname == 'p': + if k < 10 and k > 0: + path0 = reference_path + '0' + \ + str(indexes0[k]) + '/'+refname+'.h5' + + #zero_point = [12.79353, 14.32866, 6.51101] + #zero_point = np.array(zero_point) + #ndim = W.mesh().topology().dim() + #zero_point = W.tabulate_dof_coordinates().reshape((-1, ndim))[0] + + if 'cib' in mode: + barye2mmHg = 0.00750062 + else: + barye2mmHg = 1/1333.22387415 + p = Function(W) + p0 = Function(W) + + p.rename('error_p', outname) + + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(p, 'p/vector_0') + hdf0 = HDF5File(MESH['mesh'].mpi_comm(), path0, 'r') + hdf0.read(p0, 'p/vector_0') + + hdf.close() + hdf0.close() + + #p0.vector()[:] -= p0(zero_point) + pvec = p.vector().get_local() + p0vec = p0.vector().get_local() + + pvec = pvec - np.mean(pvec) + p0vec = p0vec - np.mean(p0vec) + + errvec = np.sqrt((pvec - p0vec)**2)*barye2mmHg + #errvec = np.abs((pvec - p0vec)/(p0vec)) + + p.vector()[:] = errvec + #p.vector()[:] = p0vec + + xdmf_p.write(p, te) + te = te + dt + numt = numt + 1 + +def SEQCIBH5(pathtocib, pathtocibseq, ndata): + + import csv + + times = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, + 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25] + #times = [8] + + if ndata == 1: + pathtocib += '9_mm_AoCo_phantom_rest_vox0.9mm/' + pathtocibseq += '9AoCo_rest_vox0.9/' + elif ndata == 2: + pathtocib += '9_mm_AoCo_phantom_stress_vox0.9mm/' + pathtocibseq += '9AoCo_stress_vox0.9/' + elif ndata == 3: + pathtocib += '11_mm_AoCo_phantom_rest_vox0.9mm/' + pathtocibseq += '11AoCo_rest_vox0.9/' + elif ndata == 4: + pathtocib += '11_mm_AoCo_phantom_stress_vox0.9mm/' + pathtocibseq += '11AoCo_stress_vox0.9/' + else: + raise Exception('Data file not recognize!') + + xdmf_v = XDMFFile(pathtocibseq + 'vel.xdmf') + + mesh = Mesh(pathtocib+'aorta.xml') + Evel = VectorElement('Lagrange', mesh.ufl_cell(), 1) + V = FunctionSpace(mesh, Evel) + boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1) + + VX = V.sub(0).collapse() + VY = V.sub(1).collapse() + VZ = V.sub(2).collapse() + vv_x = Function(VX) + vv_y = Function(VY) + vv_z = Function(VZ) + v = Function(V) + + v.rename('velocity', 'vel') + mx = dof_to_vertex_map(VX) + my = dof_to_vertex_map(VY) + mz = dof_to_vertex_map(VZ) + + for ti in times: + print('reading aorta '+str(ti) + '.csv') + with open(pathtocib + 'aorta'+str(ti)+'.csv', 'r') as csvfile: + mylist = [row[0] for row in csv.reader(csvfile, delimiter=';')] + + Values = np.array(mylist) + file_x = np.zeros([len(Values)]) + file_y = np.zeros([len(Values)]) + file_z = np.zeros([len(Values)]) + + for l in range(len(Values)): + row = Values[l].split(',') + file_x[l] = float(row[0]) + file_y[l] = float(row[1]) + file_z[l] = float(row[2]) + + # print(np.max(file_z)) + #print(np.where(file_z==np.max(file_z)) ) + # print(np.max(file_x[np.where(file_z==np.max(file_z))])) + # print(np.max(file_y[np.where(file_z==np.max(file_z))])) + + #S = np.load(pathtocibseq) + #ux = S['x'] + #uy = S['y'] + #uz = S['z'] + #ux = ux.transpose((0,2,1,3)) + #uy = uy.transpose((0,2,1,3)) + #uz = uz.transpose((0,2,1,3)) + + vv_x.vector()[:] = file_x[mx] + vv_y.vector()[:] = file_y[my] + vv_z.vector()[:] = file_z[mz] + + assign(v.sub(0), vv_x) + assign(v.sub(1), vv_y) + assign(v.sub(2), vv_z) + xdmf_v.write(v, ti) + # LagrangeInterpolator.interpolate(vv_x,v1_x) + # LagrangeInterpolator.interpolate(vv_y,v1_y) + # LagrangeInterpolator.interpolate(vv_z,v1_z) + + xdmf_v.close() + +def SqtoH5(BOX, MESH, Sqx, Sqy, Sqz, output_path, uname): + + xdmf_u = XDMFFile(output_path+uname+'.xdmf') + Xt = BOX['nodes'] + SqXt = BOX['seq'] + P1 = BOX['FEM'] + V = MESH['FEM'] + [Lx, Ly, Lz, Lt] = Sqx.shape + + VX = V.sub(0).collapse() + VY = V.sub(1).collapse() + VZ = V.sub(2).collapse() + + vv_x = Function(VX) + vv_y = Function(VY) + vv_z = Function(VZ) + + m = dof_to_vertex_map(P1) + + if rank == 0: + print('{SQtoH5} total number of timesteps: ' + str(Lt)) + + v = Function(V) + v.rename('velocity', uname) + dt = (Lt-1) + te = 0 + + for t in range(Lt): + + if rank == 0: + print('timestep number', t) + v1_x = Function(P1) + v1_y = Function(P1) + v1_z = Function(P1) + + values_x = np.zeros(v1_x.vector().get_local().shape) + values_y = np.zeros(v1_y.vector().get_local().shape) + values_z = np.zeros(v1_z.vector().get_local().shape) + + S0x = Sqx[:, :, :, t] + S0y = Sqy[:, :, :, t] + S0z = Sqz[:, :, :, t] + for k in range(Xt.shape[0]): + values_x[k] = S0x[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] + values_y[k] = S0y[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] + values_z[k] = S0z[SqXt[k, 0], SqXt[k, 1], SqXt[k, 2]] + + v1_x.vector()[:] = values_x + v1_y.vector()[:] = values_y + v1_z.vector()[:] = values_z + + LagrangeInterpolator.interpolate(vv_x, v1_x) + LagrangeInterpolator.interpolate(vv_y, v1_y) + LagrangeInterpolator.interpolate(vv_z, v1_z) + + assign(v.sub(0), vv_x) + assign(v.sub(1), vv_y) + assign(v.sub(2), vv_z) + + xdmf_u.write(v, te) + te = te+dt + +def ESTIMpressure(MESH, outpath, checkpoint_path, filename, options): + + #from dolfin import HDF5File + V = MESH['FEM'] + W = MESH['FEM'].sub(1).collapse() + # Checkpoints folders + unsort_indexes = os.listdir(checkpoint_path) + indexes = [int(x) for x in unsort_indexes] + indexes.sort() + + from pathlib import Path + import pickle + + # Create the Spheres + if options['Estim_Pressure']['method'] == 'spheres': + import mshr + Npts0 = options['Estim_Pressure']['spheres'][0]['Npts'] + Npts1 = options['Estim_Pressure']['spheres'][1]['Npts'] + center0 = options['Estim_Pressure']['spheres'][0]['center'] + center1 = options['Estim_Pressure']['spheres'][1]['center'] + radius0 = options['Estim_Pressure']['spheres'][0]['radius'] + radius1 = options['Estim_Pressure']['spheres'][1]['radius'] + x0 = Point(np.array(center0, dtype=float)) + x1 = Point(np.array(center1, dtype=float)) + sphere0 = mshr.Sphere(x0, radius0) + sphere1 = mshr.Sphere(x1, radius1) + mesh_s1 = mshr.generate_mesh(sphere0, Npts0) + mesh_s2 = mshr.generate_mesh(sphere1, Npts1) + + VS1 = FunctionSpace(mesh_s1, FiniteElement('P', mesh_s1.ufl_cell(), 1)) + VS2 = FunctionSpace(mesh_s2, FiniteElement('P', mesh_s2.ufl_cell(), 1)) + s1 = Function(VS1) + s2 = Function(VS2) + p = Function(W) + + one_mesh = interpolate(Constant(1), W) + LagrangeInterpolator.interpolate(s1, one_mesh) + LagrangeInterpolator.interpolate(s2, one_mesh) + vol1 = assemble(s1*dx) + vol2 = assemble(s1*dx) + + dt = options['Estim_Pressure']['dt'] + p_drop_lst = [] + time_ = [] + + for k in range(0, len(indexes)): + path = checkpoint_path + str(indexes[k]) + '/'+filename+'.h5' + hdf = HDF5File(MESH['mesh'].mpi_comm(), path, 'r') + hdf.read(p, 'p/vector_0') + hdf.close() + + LagrangeInterpolator.interpolate(s1, p) + LagrangeInterpolator.interpolate(s2, p) + + P1 = assemble(s1*dx)/vol1 + P2 = assemble(s2*dx)/vol2 + + p_drop_lst.append(P2-P1) + time_.append(k*dt) + if rank == 0: + print('Pressure drop :', P2-P1) + + # Saving the Result + write_path = Path(outpath) + write_path.mkdir(exist_ok=True) + methods = filename[2:] + data = { + 'pdrop': np.array(p_drop_lst), + 'time': np.array(time_) + } + fpath = write_path.joinpath('pdrop_' + methods + '.dat') + pickle.dump(data, fpath.open('wb')) + +def OUTLETwind(MESH, output_path, checkpoint_path, filename, bnds): + + #from dolfin import HDF5File + V = MESH['FEM'].sub(0).collapse() + W = MESH['FEM'].sub(1).collapse() + # Checkpoints folders + unsort_indexes = os.listdir(checkpoint_path) + indexes = [int(x) for x in unsort_indexes] + indexes.sort() + + QQ = {} + PP = {} + ds = Measure('ds', domain=MESH['mesh'], subdomain_data=MESH['boundaries']) + n = FacetNormal(MESH['mesh']) + ones = interpolate(Constant(1), W) + for bb in bnds: + PP[bb] = [] + QQ[bb] = [] + + for k in range(0, len(indexes)): + path_u = checkpoint_path + str(indexes[k]) + '/'+filename[0]+'.h5' + path_p = checkpoint_path + str(indexes[k]) + '/'+filename[1]+'.h5' + if indexes[k] < 10: + path_u = checkpoint_path + '0' + \ + str(indexes[k]) + '/'+filename[0]+'.h5' + path_p = checkpoint_path + '0' + \ + str(indexes[k]) + '/'+filename[1]+'.h5' + + u = Function(V) + p = Function(W) + #comm = MPI.COMM_WORLD + hdf_u = HDF5File(MESH['mesh'].mpi_comm(), path_u, 'r') + hdf_u.read(u, 'u/vector_0') + hdf_u.close() + hdf_p = HDF5File(MESH['mesh'].mpi_comm(), path_p, 'r') + hdf_p.read(p, 'p/vector_0') + hdf_p.close() + + print('Computing flows and pressures at timestep number ' + str(k)) + for bb in bnds: + if bb == 2: + QQ[bb].append(-assemble(dot(u, n)*ds(bb))) + else: + QQ[bb].append(assemble(dot(u, n)*ds(bb))) + + PP[bb].append(assemble(p*ds(bb))/1333.22387415 / + assemble(ones*ds(bb))) + + print('saving the figure at' + output_path) + from matplotlib import pyplot as plt + from matplotlib import rc + rc('text', usetex=True) + fig = plt.figure(figsize=(12, 5), dpi=100) + t = np.linspace(0, len(indexes), len(indexes)) + plt.subplot(1, 2, 2) + for bb in bnds: + plt.plot(t, QQ[bb], linewidth=1.2, linestyle='-', + label='$ Outlet: '+str(bb)+'$') + plt.xlabel(r'$frames $', fontsize=20) + plt.ylabel(r'$ Q $', fontsize=20) + plt.xlim([0, 1.05*max(t)]) + plt.legend(fontsize=14) + + plt.subplot(1, 2, 1) + for bb in bnds: + plt.plot(t, PP[bb], linewidth=1.2, linestyle='-', + label='$ Outlet: '+str(bb)+'$') + plt.xlabel(r'$frames $', fontsize=20) + plt.ylabel(r'$ P \ \ mmHg$', fontsize=20) + plt.xlim([0, 1.05*max(t)]) + plt.legend(fontsize=14) + + fig.savefig(output_path + 'Outlet_Windkessel', dpi=500) + +def ROUTINE(options): + + if 'Outlet_Wind' in options: + if options['Outlet_Wind']['apply']: + if rank == 0: + print('--- Outlet Windkessel ---') + + mesh_path = options['Outlet_Wind']['mesh_path'] + out_path = options['Outlet_Wind']['out_path'] + filename = options['Outlet_Wind']['filename'] + checkpoint = options['Outlet_Wind']['checkpoint'] + 'checkpoint/' + bnds = options['Outlet_Wind']['bnds'] + MESH = LOADmesh(mesh_path) + OUTLETwind(MESH, out_path, checkpoint, filename, bnds) + + if 'Corrector' in options: + if options['Corrector']['apply']: + if rank == 0: + print('Applying Corrector') + + MESH = LOADmesh(options['Corrector']['mesh_path']) + u_path = options['Corrector']['u_path'] + 'checkpoint/' + w_path = options['Corrector']['w_path'] + 'checkpoint/' + wname = options['Corrector']['wname'] + uname = options['Corrector']['uname'] + mode = options['Corrector']['mode'] + outpath = options['Corrector']['outpath'] + + ERRORmap(MESH, mode, outpath, u_path, + w_path, uname, wname, options) + + if 'Velocity' in options: + if options['Velocity']['apply']: + if rank == 0: + print('--- Reading Velocity ---') + + MESH = LOADmesh(options['Velocity']['mesh_path']) + filename = options['Velocity']['filename'] + checkpoint_path = options['Velocity']['checkpoint'] + 'checkpoint/' + outpath = options['Velocity']['checkpoint'] + WORKcheck(MESH, 'u', outpath, checkpoint_path, + filename, filename, options) + + if 'Vel_from_check' in options: + if options['Vel_from_check']['apply']: + if rank == 0: + print('Applying Velocity--MAP from sequence') + + [Sqx, Sqy, Sqz] = LOADsequences(options['Velocity']['pathtoseq']) + BOX = LOADmesh('cib') + AORTA = LOADmesh(pathtomesh) + SqtoH5(BOX, AORTA, Sqx, Sqy, Sqz, output_path, uname) + + if 'W_map' in options: + if options['W_map']['apply']: + if rank == 0: + print('Applying W--MAP') + MESH = LOADmesh(options['mesh_path']) + filename = options['W_map']['filename'] + outname = options['W_map']['out_name'] + mode = 'w' + WORKcheck(MESH, mode, output_path, checkpoint_path, + filename, outname, options) + + if 'GradW_map' in options: + if options['GradW_map']['apply']: + if rank == 0: + print('Applying Grad W--MAP') + MESH = LOADmesh(options['mesh_path']) + filename = options['GradW_map']['filename'] + outname = options['GradW_map']['out_name'] + mode = 'gradw' + WORKcheck(MESH, mode, output_path, checkpoint_path, + filename, outname, options) + + if 'Pressure' in options: + if options['Pressure']['apply']: + if rank == 0: + print('Applying Pressure-MAP') + + MESH = LOADmesh(options['mesh_path']) + filename = options['Pressure']['filename'] + outname = options['Pressure']['out_name'] + mode = 'p' + WORKcheck(MESH, mode, output_path, checkpoint_path, + filename, outname, options) + + if 'Error_P' in options: + if options['Error_P']['apply']: + if rank == 0: + print('Applying L2 error to Pressure') + + MESH = LOADmesh(options['mesh_path']) + refname = options['Error_P']['refname'] + reference_path = options['Error_P']['refpath'] + filename = options['Error_P']['filename'] + outname = options['Error_P']['out_name'] + mode = 'error_p' + ERRORmap(MESH, mode, output_path, reference_path, + checkpoint_path, refname, filename, outname, options) + + if 'Estim_Pressure' in options: + if options['Estim_Pressure']['apply']: + if rank == 0: + print('Applying Pressure Estimator') + + MESH = LOADmesh(options['mesh_path']) + filename = options['Estim_Pressure']['filename'] + outpath = options['Estim_Pressure']['outpath'] + ESTIMpressure(MESH, outpath, checkpoint_path, filename, options) + + if 'Error-curves' in options: + if options['Error-curves']['apply']: + if rank == 0: + print('--- Error curves ---' ) + + MESH = LOADmesh(options['Error-curves']['meshpath']) + ref_check = options['Error-curves']['ref_check'] + 'checkpoint/' + com_check = options['Error-curves']['com_check'] + 'checkpoint/' + outpath = options['Error-curves']['outpath'] + refname = options['Error-curves']['ref_name'] + comname = options['Error-curves']['com_name'] + ERRORmap(MESH, 'curves', outpath, ref_check, com_check, + refname, comname, options) + + if 'Temporal-Average' in options: + if options['Temporal-Average']['apply']: + if rank == 0: + print('--- Temporal Average ---') + + MESH = LOADmesh(options['Temporal-Average']['meshpath']) + ref_check = options['Temporal-Average']['original_check'] + 'checkpoint/' + out_check = options['Temporal-Average']['out_check'] + READcheckpoint(MESH,'average', out_check,ref_check,'u','u',options) + + if 'Temporal-Interpolation' in options: + if options['Temporal-Interpolation']['apply']: + if rank == 0: + print('--- Temporal Interpolation ---') + + MESH = LOADmesh(options['Temporal-Interpolation']['meshpath']) + ref_check = options['Temporal-Interpolation']['original_check'] + 'checkpoint/' + out_check = options['Temporal-Interpolation']['out_check'] + READcheckpoint(MESH,'interpolation', out_check,ref_check,'u','u',options) + + + + +if __name__ == '__main__': + + comm = MPI.COMM_WORLD + size = comm.Get_size() + rank = comm.Get_rank() + + if len(sys.argv) > 1: + if os.path.exists(sys.argv[1]): + inputfile = sys.argv[1] + if rank == 0: + print('Found input file ' + inputfile) + else: + raise Exception('Command line arg given but input file does not exist:' + ' {}'.format(sys.argv[1])) + else: + raise Exception('An input file is required as argument!') + + options = inout.read_parameters(inputfile) + ROUTINE(options) diff --git a/codes/SENSE.py b/codes/SENSE.py new file mode 100644 index 0000000..83b82f2 --- /dev/null +++ b/codes/SENSE.py @@ -0,0 +1,115 @@ +import numpy as np +from numpy import linalg as LA +import sys +from mpi4py import MPI +comm = MPI.COMM_WORLD +size = comm.Get_size() +rank = comm.Get_rank() + +# SENSE: Simulation of SENSitive Encoding algorithm proposed by K. Pruessmann, et. al. in: +# "SENSE: Sensitivity Enconding for Fast MRI" Mag. Res. in Medicine 42. (1999) +# written by Jeremias Garay (j.e.garay.labra@rug.nl) + +def Sensitivity_Map(shape): + + [Nx,Ny,Nz] = shape + [X,Y,Z] = np.meshgrid(np.linspace(0,Ny,Ny),np.linspace(0,Nx,Nx),np.linspace(0,Nz,Nz)) + Xsense1 = (X/(Nx*2)-1)**2 + Xsense2 = ((Nx-X)/(Nx*2)-1)**2 + S_MAPS = [np.fft.fftshift(Xsense1),np.fft.fftshift(Xsense2)] + + return S_MAPS + +def SENSE_recon(S1,M1,S2,M2): + + [Nx,Ny,Nz,Nt] = M1.shape + M = np.zeros([Nx,int(2*Ny),Nz,Nt],dtype=complex) + sm1 = np.fft.fftshift(S1)[:,:,0] + sm2 = np.fft.fftshift(S2)[:,:,0] + + for j in range(Ny): + for k in range(Nx): + l1 = M1[k,j,:,:]; a1 = sm1[k,j]; a2 = sm1[k,j+Ny] + l2 = M2[k,j,:,:]; b1 = sm2[k,j]; b2 = sm2[k,j+Ny] + B = (l1*b1 - l2*a1)/(a2*b1 - b2*a1) + A = (l1*b2 - l2*a2)/(a1*b2 - a2*b1) + M[k,j,:,:] = A + M[k,j+Ny,:,:] = B + + + return M + +def SENSE_recon2(S1,M1,S2,M2): + # With matrices as in the original paper! + + [Nx,Ny,Nz,Nt] = M1.shape + M = np.zeros([Nx,int(2*Ny),Nz,Nt],dtype=complex) + sm1 = np.fft.fftshift(S1)[:,:,0] + sm2 = np.fft.fftshift(S2)[:,:,0] + sigma2 = 0.049**2 + sigma2 = 1 + Psi = np.diagflat(np.array([sigma2,sigma2])) # Error matrix Psi + Psi_inv = np.linalg.inv(Psi) + + for j in range(Ny): + for k in range(Nx): + l1 = M1[k,j,:,:]; a1 = sm1[k,j]; a2 = sm1[k,j+Ny] + l2 = M2[k,j,:,:]; b1 = sm2[k,j]; b2 = sm2[k,j+Ny] + S = np.array([[a1,a2],[b1,b2]]) + U = np.linalg.inv((np.transpose(S)*Psi_inv*S))*np.transpose(S)*Psi_inv + a = np.array([l1,l2]) + a_resized = np.resize(a,(2,Nz*Nt)) + v_resized = np.dot(U,a_resized) + v = np.resize(v_resized,(2,Nz,Nt)) + M[k,j,:,:] = v[0,:,:] + M[k,j+Ny,:,:] = v[1,:,:] + + + return M + +def SENSE_METHOD(Seq,R): + ''' + Args: + ITOT: a numpy matrix with the full sampled (3D or 4D) dynamical data + R: the acceleration factor + ''' + + [row,col,dep,numt2] = Seq.shape + Seq_red = {} + SenseMAP = {} + [SenseMAP[0],SenseMAP[1]] = Sensitivity_Map([row,col,dep]) + + col2 = int(np.ceil(col/2)) + + for rs in range(R): + Seq_red[rs] = np.zeros([row,col2,dep,numt2],dtype=complex) + for t in range(numt2): + Kdata_0 = np.fft.fftn(Seq[:,:,:,t]) + Kdata_0 = Kdata_0*SenseMAP[rs] + Kdata_0 = Kdata_0[:,0::R,:] + Seq_red[rs][:,:,:,t] = np.fft.ifftn(Kdata_0) + + Seq_recon = SENSE_recon2(SenseMAP[0],Seq_red[0],SenseMAP[1],Seq_red[1]) + + return Seq_recon + +def undersampling(Mx,My,Mz,options): + + R = options['SENSE']['R'] + + for r in R: + if rank==0: + print('Using Acceleration Factor R = ' + str(r)) + print('applying into x component') + Mx_s = SENSE_METHOD(Mx,r) + if rank==0: + print('applying into y component') + My_s = SENSE_METHOD(My,r) + if rank==0: + print('applying into z component') + Mz_s = SENSE_METHOD(Mz,r) + + return [Mx_s,My_s,Mz_s] + + + diff --git a/codes/__pycache__/CS.cpython-36.pyc b/codes/__pycache__/CS.cpython-36.pyc new file mode 100644 index 0000000..c0eede0 Binary files /dev/null and b/codes/__pycache__/CS.cpython-36.pyc differ diff --git a/codes/__pycache__/Graphics.cpython-36.pyc b/codes/__pycache__/Graphics.cpython-36.pyc new file mode 100644 index 0000000..d0cf857 Binary files /dev/null and b/codes/__pycache__/Graphics.cpython-36.pyc differ diff --git a/codes/__pycache__/MRI.cpython-36.pyc b/codes/__pycache__/MRI.cpython-36.pyc new file mode 100644 index 0000000..eaa8bcc Binary files /dev/null and b/codes/__pycache__/MRI.cpython-36.pyc differ diff --git a/codes/__pycache__/SENSE.cpython-36.pyc b/codes/__pycache__/SENSE.cpython-36.pyc new file mode 100644 index 0000000..ff25320 Binary files /dev/null and b/codes/__pycache__/SENSE.cpython-36.pyc differ diff --git a/codes/__pycache__/ktBLAST.cpython-36.pyc b/codes/__pycache__/ktBLAST.cpython-36.pyc new file mode 100644 index 0000000..7cba033 Binary files /dev/null and b/codes/__pycache__/ktBLAST.cpython-36.pyc differ diff --git a/codes/cfd_mono.py b/codes/cfd_mono.py new file mode 100644 index 0000000..b6d17d4 --- /dev/null +++ b/codes/cfd_mono.py @@ -0,0 +1,327 @@ +from dolfin import * +import matplotlib.pyplot as plt +import numpy as np +import dolfin +from common import inout +from mpi4py import MPI +import sys +import os +# +# NAVIER STOKES PROBLEM IN THE AORTA with a MONOLITHIC SOLVER +# THIS SCRIPT INCLUDE THE 0-WINDKESSEL BOUNDARY CONDITION +# +# Written by Jeremias Garay L: j.e.garay.labra@rug.nl +# + +if '2017' in dolfin.__version__: + class MPI(MPI): + comm_world = mpi_comm_world() + + set_log_active(False) +else: + parameters["std_out_all_processes"] = False + + +def BACH(): + import os + tempof = 0.7 + semicorchea = 0.1*tempof + corchea = 0.2*tempof + negra = 0.4*tempof + blanca = 0.8*tempof + LA3 = 220 + SIb = 233.08 + SI = 246.94 + DO = 261 + REb = 277.18 + RE = 293.66 + FA = 349.23 + MIb = 311.13 + MI = 329.63 + SOL = 391 + LA = 440 + + os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (negra, DO)) + os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (negra, LA)) + os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (negra, SOL)) + os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (negra, FA)) + os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (semicorchea, MI)) + os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (semicorchea, FA)) + os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (corchea, SOL)) + os.system('play -V0 --no-show-progress --null --channels 1 synth %s sine %f' % (blanca, FA)) + +def flux(u,n,ds,idd): + Q = 0 + for k in idd: + Q += assemble(dot(u,n)*ds(k)) + return Q + +def save_outlet_data(options,Qin,Qout3,Qout4,Qout5,Qout6,Pin,Pout3,Pout4,Pout5,Pout6,dt): + # saving the fluxes + np.savetxt(options['outlets_path'] + 'fluxes/' + 'q_in_dt' + str(dt) + '.txt', Qin) + np.savetxt(options['outlets_path'] + 'fluxes/' + 'q3_dt' + str(dt) + '.txt', Qout3) + np.savetxt(options['outlets_path'] + 'fluxes/' + 'q4_dt' + str(dt) + '.txt', Qout4) + np.savetxt(options['outlets_path'] + 'fluxes/' + 'q5_dt' + str(dt) + '.txt', Qout5) + np.savetxt(options['outlets_path'] + 'fluxes/' + 'q6_dt' + str(dt) + '.txt', Qout6) + # saving the pressures + np.savetxt(options['outlets_path'] + 'pressures/' + 'p_in_dt' + str(dt) + '.txt', Pin) + np.savetxt(options['outlets_path'] + 'pressures/' + 'p3_dt' + str(dt) + '.txt', Pout3) + np.savetxt(options['outlets_path'] + 'pressures/' + 'p4_dt' + str(dt) + '.txt', Pout4) + np.savetxt(options['outlets_path'] + 'pressures/' + 'p5_dt' + str(dt) + '.txt', Pout5) + np.savetxt(options['outlets_path'] + 'pressures/' + 'p6_dt' + str(dt) + '.txt', Pout6) + +def windkessel_update(u0,n,ds,fluxes,press,pii,pii0,windkessel): + # Updating the time dependent windkessel parameters + if windkessel['C']: + for nk in windkessel['id']: + k = str(nk) + fluxes[k].assign(assemble(dot(u0,n)*ds(nk))) + pii0[k].assign(pii[k]) + pii[k].assign(windkessel['alpha'][k]*pii0[k] + windkessel['beta'][k]*fluxes[k]) + press[k].assign(windkessel['gamma'][k]*fluxes[k] + windkessel['alpha'][k]*pii0[k]) + + else: + for nk in windkessel['id']: + k = str(nk) + fluxes[k] = assemble(dot(u0,n)*ds(nk)) + press[k].assign( windkessel['R_p'][k]*assemble(dot(u0,n)*ds(nk))) + +def solv_NavierStokes(options): + + # Assign physical parameters + rho = Constant(options['density']) + mu = Constant(options['dynamic_viscosity']) + otheta = Constant(1-options['param']['theta']) + theta = Constant(options['param']['theta']) + Tf = options['Tf'] + dt = options['dt'] + dt_w = options['dt_write'] + Qin = np.zeros([int(Tf/dt)]) + Qout3 = np.zeros([int(Tf/dt)]) + Qout4 = np.zeros([int(Tf/dt)]) + Qout5 = np.zeros([int(Tf/dt)]) + Qout6 = np.zeros([int(Tf/dt)]) + Pin = np.zeros([int(Tf/dt)]) + Pout3 = np.zeros([int(Tf/dt)]) + Pout4 = np.zeros([int(Tf/dt)]) + Pout5 = np.zeros([int(Tf/dt)]) + Pout6 = np.zeros([int(Tf/dt)]) + barye2mmHg = 1/1333.22387415 + + # CREATING THE FILES + xdmf_u = XDMFFile(options['savepath']+'u.xdmf') + xdmf_p = XDMFFile(options['savepath']+'p.xdmf') + xdmf_u.parameters['rewrite_function_mesh'] = False + xdmf_p.parameters['rewrite_function_mesh'] = False + # LOADING THE MESH + #mesh = Mesh(options['mesh_path']) + #boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1) + #boundaries = MarkBoundaries(boundaries) + + mesh = Mesh() + hdf = HDF5File(mesh.mpi_comm(), options['mesh_path'] , 'r') + hdf.read(mesh, '/mesh', False) + boundaries = MeshFunction('size_t', mesh , mesh.topology().dim() - 1) + hdf.read(boundaries, '/boundaries') + + # To save the boundaries information + #path2 = '/home/p283370/Desktop/marked_fine/boundaries.xdmf' + #XDMFFile(path2).write(boundaries) + # DEFINE FUNCTION SPACES + V = VectorElement('Lagrange', mesh.ufl_cell(), options['param']['Nvel']) + Q = FiniteElement('Lagrange', mesh.ufl_cell(), options['param']['Npress']) + TH = V * Q + W = FunctionSpace(mesh, TH) + # No-slip boundary condition for velocity on walls + noslip = Constant((0, 0, 0)) + inflow = Expression(('0','0','(-0.5*U*fabs(sin(w*t)) - 0.5*U*sin(w*t))*(t<0.7)'), degree=3,t=0,U=options['param']['U'], w=2*DOLFIN_PI/options['param']['period']) + bc_inflow = DirichletBC(W.sub(0), inflow, boundaries, 2) + bc_walls = DirichletBC(W.sub(0), noslip, boundaries, 1) + # COLLECTING THE BOUNDARY CONDITIONS + BCS = [bc_inflow,bc_walls] + u, p = TrialFunctions(W) + v, q = TestFunctions(W) + w = Function(W) + n = FacetNormal(mesh) + ds = Measure("ds", subdomain_data=boundaries) + h = CellDiameter(mesh) + f = Constant((0,0,0)) + u0, p0 = w.split() + + u_ = theta*u + (1 - theta)*u0 + theta_p = theta + p_ = theta_p*p + (1 - theta_p)*p0 + # The variational formulation + # Mass matrix + F = ( + (rho/dt)*dot(u - u0, v)*dx + + mu*inner(grad(u_), grad(v))*dx + - p_*div(v)*dx + q*div(u)*dx + + rho*dot(grad(u_)*u0, v)*dx + ) + + if options['stab']['temam']: + F += 0.5*rho*div(u0)*dot(u_, v)*dx + + if options['stab']['pspg']: + eps = Constant(options['stab']['epsilon']) + F += eps/mu*h**2*inner(grad(p_), grad(q))*dx + + if options['stab']['backflow']: + def abs_n(x): + return 0.5*(x - abs(x)) + + for nk in options['stab']['back_id']: + if rank==0: + print('adding backflow stabilization in border number:' + str(nk)) + F -= 0.5*rho*abs_n(dot(u0, n))*dot(u_, v)*ds(nk) + + a = lhs(F) + L = rhs(F) + + # Initialization of Windkessel Boundaries + if options['windkessel']['id']: + windkessel = options['windkessel'] + # Coeficients + fluxes = {str(k):[] for k in windkessel['id']} + pii0 = {str(k):[] for k in windkessel['id']} + pii = {str(k):[] for k in windkessel['id']} + press = {str(k):[] for k in windkessel['id']} + + if windkessel['C']: + for nk in windkessel['id']: + k = str(nk) + if rank==0: + print('Capacitance of windkessel model is: ', windkessel['C'][k]) + + # computing coeficients + windkessel['alpha'][k] = windkessel['R_d'][k]*windkessel['C'][k]/(windkessel['R_d'][k]*windkessel['C'][k] + options['dt']) + windkessel['beta'][k] = windkessel['R_d'][k]*(1-windkessel['alpha'][k]) + windkessel['gamma'][k] = windkessel['R_p'][k] + windkessel['beta'][k] + if rank==0: + print('Using 0-Windkessel complete at outlet number: ' + str(k)) + # setting initial values for flux, distal pressure and proximal pressure + fluxes[k] = Constant(0) + pii0[k] = Constant(47/barye2mmHg) + pii[k] = Constant(windkessel['alpha'][k]*pii0[k] + windkessel['beta'][k]*fluxes[k]) + press[k] = Constant(windkessel['gamma'][k]*fluxes[k] + windkessel['alpha'][k]*pii0[k]) + # Adding to RHS + L = L - dt*press[k]*dot(v,n)*ds(nk) + else: + for nk in windkessel['id']: + k = str(nk) + if rank==0: + print('Using 0-Windkessel reduced at outlet number: ' + str(nk)) + fluxes[k] = Constant(0) + press[k] = Constant(windkessel['R_p'][k]*0) + # Adding to RHS + L = L - dt*press[k]*dot(v,n)*ds(nk) + + + + # The static part of the matrix + A = assemble(a) + #b = assemble(L) + [bc.apply(A) for bc in BCS] + #[bc.apply(b) for bc in BCS] + #solv = LUSolver() + #solv.set_operator(A) + #solv.parameters['linear_solver'] = 'mumps' + #solv.parameters['reuse_factorization'] = True + u, p = w.split() + u.rename('velocity', 'u') + p.rename('pressure', 'p') + ind = 0 + t = dt + + ones = interpolate(Constant(1),W.sub(1).collapse()) + A2 = assemble(ones*ds(2)) + A3 = assemble(ones*ds(3)) + A4 = assemble(ones*ds(4)) + A5 = assemble(ones*ds(5)) + A6 = assemble(ones*ds(6)) + + checkcicle = int(options['checkpoint_dt']/options['dt']) + writecicle = int(options['checkpoint_dt']/options['dt_write']) + + while t<=Tf+dt: + + if options['windkessel']['id']: + windkessel_update(u,n,ds,fluxes,press,pii,pii0,windkessel) + # To solve + assemble(a, tensor=A) + b = assemble(L) + [bc.apply(A, b) for bc in BCS] + solve(A, w.vector(), b) + + Qin[ind] = flux(u,n,ds,[2]) + Qout3[ind] = flux(u,n,ds,[3]) + Qout4[ind] = flux(u,n,ds,[4]) + Qout5[ind] = flux(u,n,ds,[5]) + Qout6[ind] = flux(u,n,ds,[6]) + Pin[ind] = barye2mmHg*assemble(p*ds(2))/A2 + Pout3[ind] = barye2mmHg*assemble(p*ds(3))/A3 + Pout4[ind] = barye2mmHg*assemble(p*ds(4))/A4 + Pout5[ind] = barye2mmHg*assemble(p*ds(5))/A5 + Pout6[ind] = barye2mmHg*assemble(p*ds(6))/A6 + + if rank==0: + print('t = ',t) + # print('|u|:', norm(u0)) + # print('|p|:', norm(p0)) + # print('div(u):', assemble(div(u0)*dx)) + print('Dp = ',np.round(Pin[ind]-Pout3[ind],3),'mmHg') + + ind += 1 + if options['write_xdmf']: + if np.mod(ind,writecicle)<0.1 or ind==1: + xdmf_u.write(u, t) + xdmf_p.write(p, t) + + if np.mod(ind,checkcicle)<0.1 or ind==1: + if options['write_checkpoint']: + checkpath = options['savepath'] +'checkpoint/{i}/'.format(i=ind) + comm = u.function_space().mesh().mpi_comm() + inout.write_HDF5_data(comm, checkpath + '/u.h5', u, '/u', t=t) + inout.write_HDF5_data(comm, checkpath + '/p.h5', p, '/p', t=t) + + inflow.t = t + t += dt + + + # saving the data at outlets: fluxes and pressures + if options['write_outlets']: + if rank==0: + save_outlet_data(options,Qin,Qout3,Qout4,Qout5,Qout6,Pin,Pout3,Pout4,Pout5,Pout6,options['dt']) + + +if __name__ == '__main__': + + comm = MPI.COMM_WORLD + size = comm.Get_size() + rank = comm.Get_rank() + + if len(sys.argv) > 1: + if os.path.exists(sys.argv[1]): + inputfile = sys.argv[1] + if rank==0: + print('Found input file ' + inputfile) + else: + raise Exception('Command line arg given but input file does not exist:' + ' {}'.format(sys.argv[1])) + else: + raise Exception('An input file is required as argument!') + + if 'Zion' in os.popen('hostname').read(): + user = 'yeye' + np.set_printoptions(threshold=5) + + if 'fwn-bborg-5-166' in os.popen('hostname').read(): + user = 'p283370' + + if rank==0: + print('Welcome user {uss}'.format(uss=user)) + + options = inout.read_parameters(inputfile) + solv_NavierStokes(options) + diff --git a/codes/ktBLAST.py b/codes/ktBLAST.py new file mode 100644 index 0000000..7dadb9c --- /dev/null +++ b/codes/ktBLAST.py @@ -0,0 +1,870 @@ +import numpy as np +import scipy as sc +from scipy import signal +from mpi4py import MPI +comm = MPI.COMM_WORLD +size = comm.Get_size() +rank = comm.Get_rank() + + + +# kt-BLAST (NO DC TERM) method for reconstruction of undersampled MRI image based on +# l2 minimization. + +def EveryAliased3D2(i,j,k,PP,Nx,Ny,Nz,BB,R): + + ivec = [i,j,k] + Nvec = [Nx,Ny,Nz] + [ktot,ltot] = PP.shape + Ptot = np.zeros([ktot**ltot,ltot]) + PP2 = np.zeros([ktot**ltot,ltot]) + tt = -1 + + for kk in range(Ptot.shape[0]): + nn = int(np.mod(kk,3)) + mm = int(np.mod(np.floor(kk/3),3)) + if np.mod(kk,9)==0: + tt+=1 + + Ptot[kk,0] = PP[tt,0] + ivec[0] + Ptot[kk,1] = PP[mm,1] + ivec[1] + Ptot[kk,2] = PP[nn,2] + ivec[2] + + for kk in range(Ptot.shape[0]): + for ll in range(Ptot.shape[1]): + if Ptot[kk,ll]<0: + Ptot[kk,ll] = Ptot[kk,ll] + Nvec[ll] + if Ptot[kk,ll]>=Nvec[ll]: + Ptot[kk,ll] = Ptot[kk,ll] - Nvec[ll] + + + CC = np.zeros([3,Ptot.shape[0]+1]) + YY = np.array([ [i] , [j], [k] ]) + CC[0,0] = i + CC[1,0] = j + CC[2,0] = k + psel = 0 + + + for l in range(1,Ptot.shape[0]+1): + CC[0,l] = int(Ptot[l-1,0]) + CC[1,l] = int(Ptot[l-1,1]) + CC[2,l] = int(Ptot[l-1,2]) + + + + if CC[0,l]==YY[0,psel] and CC[1,l]==YY[1,psel] and CC[2,l]==YY[2,psel] and BB[int(CC[1,l]),int(CC[2,l]),int(CC[0,l])]!=0: + pass + else: + War = False + for ww in range(psel): + if CC[0,l]==YY[0,ww] and CC[1,l]==YY[1,ww] and CC[2,l]==YY[2,ww] and BB[int(CC[1,l]),int(CC[2,l]),int(CC[0,l])]!=0: + War = True + + if not War: + psel += 1 + CCC = np.array([ [CC[0,l] ] , [CC[1,l]] , [CC[2,l]]]) + YY = np.concatenate( ( YY, CCC ) ,axis=1 ) + + return YY.astype(int) + +def EveryAliased3D(i,j,k,DP,Nx,Ny,Nz,BB,R,SPREAD=None): + + ivec = [i,j,k] + Nvec = [Nx,Ny,Nz] + [ktot,ltot] = DP.shape + DPN = np.zeros([ktot,ltot]) + + if SPREAD is not None: # WITH SPREAD FUNCTIONS FORMALISM + Maux = np.zeros([Ny,Nz,Nx]) + Maux[j,k,i] = 1 + SP2 = SPREAD[::-1,::-1,::-1] + MS = R*sc.signal.convolve(Maux,SP2, mode='same') + ms = np.abs(MS) + Ims = 1*(ms>np.max(ms)*0.405) + Pas = np.where(Ims==1) + PP = np.array(Pas[:]) + PEA = PP[::-1,:] + + for ll in range(PEA.shape[1]): + if PEA[0,ll]>=Nx: + PEA[0,ll] = PEA[0,ll] - Nx + if PEA[1,ll]>=Ny: + PEA[1,ll] = PEA[1,ll] - Ny + if PEA[2,ll]>=Nz: + PEA[2,ll] = PEA[2,ll] - Nz + + + Ntot = PEA.shape[1] + ind = 0 + PEAnew = PEA + + + for ll in range(Ntot): + if BB[PEA[1,ll],PEA[2,ll],PEA[0,ll]]!=0: + PEAnew = np.delete(PEAnew,(ll-ind),axis=1) + ind +=1 + + + return PEA + else: + + for kk in range(DPN.shape[0]): + for l in range(DPN.shape[1]): + DPN[kk,l] = DP[kk,l] + ivec[l] + if DPN[kk,l]<0: + DPN[kk,l] = DPN[kk,l] + Nvec[l] + if DPN[kk,l]>=Nvec[l]: + DPN[kk,l] = DPN[kk,l] - Nvec[l] + + + + CC = np.zeros([3,ktot+1]) + YY = np.array([ [i] , [j], [k] ]) + CC[0,0] = i + CC[1,0] = j + CC[2,0] = k + + + for l in range(1,ktot+1): + CC[0,l] = DPN[l-1,0] + CC[1,l] = DPN[l-1,1] + CC[2,l] = DPN[l-1,2] + + if CC[0,l]!=CC[0,l-1] and CC[1,l]!=CC[1,l-1] and CC[2,l]!=CC[2,l-1] and BB[int(CC[1,l]),int(CC[2,l]),int(CC[0,l])]==0: + CCC = np.array([ [CC[0,l] ] , [CC[1,l]] , [CC[2,l]]]) + YY = np.concatenate( ( YY, CCC ) ,axis=1 ) + + return YY.astype(int) + +def EveryAliased(i,j,DP,Nx,Ny,BB,R,mode): + + if mode==1: # USING GEOMETRICAL ASSUMPTIONS + ivec = [i,j] + Nvec = [Nx,Ny] + DPN = 0*DP + [ktot,ltot] = DP.shape + + for k in range(ktot): + for l in range(ltot): + + DPN[k,l] = DP[k,l] + ivec[l] + + if DPN[k,l]<0: + #DPN[k,l] = ivec[l] + DPN[k,l] = DPN[k,l] + Nvec[l] + + if DPN[k,l]>=Nvec[l]: + #DPN[k,l] = ivec[l] + DPN[k,l] = DPN[k,l] - Nvec[l] + + CC = np.zeros([2,ktot+1]) + YY = np.array([ [i] , [j] ]) + CC[0,0] = i + CC[1,0] = j + + for l in range(1,ktot+1): + CC[0,l] = DPN[l-1,0] + CC[1,l] = DPN[l-1,1] + if CC[0,l]!=CC[0,l-1] and CC[1,l]!=CC[1,l-1] and BB[int(CC[1,l]),int(CC[0,l])]==0: + CCC = np.array([ [CC[0,l] ] , [CC[1,l]] ]) + YY = np.concatenate( ( YY, CCC ) ,axis=1 ) + + return YY.astype(int) + + + + if mode=='spread': # WITH SPREAD FUNCTIONS FORMALISM + Maux = np.zeros([row,numt2]) + Maux[l,k] = 1 + MS = R*ConvolSP(Maux,SPREAD) + ms = np.abs(MS) + Ims = 1*(ms>np.max(ms)*0.405) + Pas = np.where(Ims==1) + PP = np.array(Pas[:]) + return PP[::-1,:] + +def GetSymmetric(M): + [row,numt2] = M.shape + S = np.zeros(M.shape,dtype=complex) + aux = np.zeros([1,row]) + nmid = 0.5*(numt2+1) + for k in range(int(nmid)): + aux = 0.5*( M[:,k] + M[:,numt2-k-1] ) + S[:,k] = aux + S[:,numt2-k-1] = aux + + return S + +def UNDER(A,mode,R,k): + + start = np.mod(k,R) + I1B = np.zeros(A.shape,dtype=complex) + + # Not quite efficient ! better to work with vectors + + if mode=='ky': + for k in range(start,A.shape[0],R): + I1B[k,:,:] = A[k,:,:] + + if mode=='kxky': + for k in range(start,A.shape[0],R): + for l in range(start,A.shape[1],R): + I1B[k,l,:] = A[k,l,:] + + if mode=='kxkykz': + for k in range(start,A.shape[0],R): + for l in range(start,A.shape[2],R): + for r in range(start,A.shape[1],R): + I1B[k,r,l] = A[k,r,l] + + + return I1B + +def FilteringHigh(M,fac): + + if M.ndim==2: + + [row,col] = M.shape + inx = np.linspace(0,col-1,col) + MF = np.zeros(M.shape,dtype=complex) + + for k in range(row): + vecfou = np.fft.fft(M[k,:]) + window = signal.tukey(2*col,fac) + vecfou2 = vecfou*window[col:2*col] + MF[k,:] = np.fft.ifft(vecfou2) + + return MF + + + if M.ndim==3: + [row,col,dep] = M.shape + MF = np.zeros(M.shape,dtype=complex) + inx = np.linspace(0,col-1,col) + for l in range(dep): + for k in range(row): + vecfou = np.fft.fft(M[k,:,l]) + window = signal.tukey(2*col,fac) + vecfou2 = vecfou*window[col:2*col] + MF[k,:,l] = np.fft.ifft(vecfou2) + + + return MF + +def InterpolateM(M,numt2): + + if M.ndim==2: + + [row,numt] = M.shape + MNew = np.zeros([row,numt2],dtype=complex) + xdat = np.linspace(0,numt,numt) + xdat2 = np.linspace(0,numt,numt2) + nstar = int(0.5*(numt2-numt)) + + for t in range(nstar,nstar+numt): + MNew[:,t] = M[:,t-nstar] + + for l in range(row): + ydat = M[l,:] + fdat = sc.interpolate.interp1d(xdat,ydat,kind='cubic') + MNew[l,1:nstar] = fdat(xdat2)[1:nstar] + MNew[l,nstar+numt:numt2] = fdat(xdat2)[nstar+numt:numt2] + + + if M.ndim==3: + [row,col,numt] = M.shape + MNew = np.zeros([row,col,numt2],dtype=complex) + xdat = np.linspace(0,numt,numt) + xdat2 = np.linspace(0,numt,numt2) + nstar = int(0.5*(numt2-numt)) + for t in range(nstar,nstar+numt): + MNew[:,:,t] = M[:,:,t-nstar] + + for c in range(col): + for l in range(row): + ydat = M[l,c,:] + fdat = sc.interpolate.interp1d(xdat,ydat,kind='cubic') + MNew[l,c,1:nstar] = fdat(xdat2)[1:nstar] + MNew[l,c,nstar+numt:numt2] = fdat(xdat2)[nstar+numt:numt2] + + + return MNew + +def KTT(M,scol): + + #Maux = M[:,scol,:] + #Maux = np.fft.ifftshift(Maux,axes=0) + #Maux = np.fft.ifft(Maux,axis=0) + #Maux = np.fft.ifft(Maux,axis=1) + #Maux = np.fft.fftshift(Maux,axes=1) + + # TAO STYLE + Maux = np.zeros(M.shape,dtype=complex) + + for k in range(M.shape[2]): + Maux[:,:,k] = np.fft.ifftshift(M[:,:,k]) + Maux[:,:,k] = np.fft.ifft2(Maux[:,:,k]) + Maux = Maux[:,scol,:] + Maux = np.fft.ifft(Maux,axis=1) + Maux = np.fft.fftshift(Maux,axes=1) + return Maux + +def IKTT(M): + + #Maux = np.fft.ifftshift(M,axes=1) + #Maux = np.fft.ifft(Maux,axis=1) + #Maux = np.fft.fft(Maux,axis=0) + #Maux = np.fft.fftshift(Maux,axes=0) + + # TAO STYLE + Maux = np.fft.ifftshift(M,axes=1) + Maux = np.fft.fft(Maux,axis=1) + + return Maux + +def KTT3D(M,sdep): + + Maux = np.zeros(M.shape,dtype=complex) + for k in range(M.shape[3]): + Maux[:,:,:,k] = np.fft.ifftshift(M[:,:,:,k]) + Maux[:,:,:,k] = np.fft.ifftn(Maux[:,:,:,k]) + Maux = Maux[:,:,sdep,:] + Maux = np.fft.ifft(Maux,axis=2) + Maux = np.fft.fftshift(Maux,axes=2) + return Maux + +def IKTT3D(M): + + Maux = np.fft.ifftshift(M,axes=2) + Maux = np.fft.fft(Maux,axis=2) + + return Maux + +def get_points4D(row,col,dep,numt2,R,mode): + + bb = np.ceil(row/R) + cc = np.ceil(col/R) + aa = np.ceil(numt2/R) + points = R+1 + kmid = int(R/2) + + if mode=='kxky': + PC = [np.ceil(numt2/2),np.ceil(row/2),np.ceil(col/2)] + PP = np.zeros([points,3]) + DP = np.zeros([points-1,3]) + for k in range(points): + + PP[k,0] = numt2-aa*(k) + 1 + PP[k,1] = bb*(k) + PP[k,2] = cc*(k) + + if PP[k,0]>=numt2: + PP[k,0] -= 1 + if PP[k,0]<0: + PP[k,0] += 1 + + if PP[k,1]>=row: + PP[k,1] -= 1 + if PP[k,1]<0: + PP[k,1] += 1 + + if PP[k,2]>=col: + PP[k,2] -= 1 + if PP[k,2]<0: + PP[k,2] += 1 + + if kkmid: + DP[k-1,0] = PP[k,0] - PC[0] + DP[k-1,1] = PP[k,1] - PC[1] + DP[k-1,2] = PP[k,2] - PC[2] + + kmax = int((PP[kmid,0] + PP[kmid-1,0])/2 ) + kmin = int((PP[kmid,0] + PP[kmid+1,0])/2 ) + cmax = int((PP[kmid,1] + PP[kmid-1,1])/2 ) + cmin = int((PP[kmid,1] + PP[kmid+1,1])/2 ) + + + #DP2 = np.zeros([DP.shape[0]**DP.shape[1],DP.shape[1]]) + #DP2[0,0] = DP[0,0]; DP2[0,1] = DP[0,1] ; DP2[0,2] = DP[0,2] + #DP2[1,0] = DP[0,0]; DP2[1,1] = DP[0,1] ; DP2[1,2] = DP[1,2] + #DP2[2,0] = DP[0,0]; DP2[2,1] = DP[1,1] ; DP2[2,2] = DP[0,2] + #DP2[3,0] = DP[0,0]; DP2[3,1] = DP[1,1] ; DP2[3,2] = DP[1,2] + #DP2[4,0] = DP[1,0]; DP2[4,1] = DP[0,1] ; DP2[4,2] = DP[0,2] + #DP2[5,0] = DP[1,0]; DP2[5,1] = DP[0,1] ; DP2[5,2] = DP[1,2] + #DP2[6,0] = DP[1,0]; DP2[6,1] = DP[1,1] ; DP2[6,2] = DP[0,2] + #P2[7,0] = DP[1,0]; DP2[7,1] = DP[1,1] ; DP2[7,2] = DP[1,2] + + + return [kmin,kmax,PP,DP] + + if mode=='ky': + PC = [np.ceil(numt2/2),np.ceil(row/2)] + PP = np.zeros([points,2]) + DP = np.zeros([points-1,2]) + for k in range(points): + PP[k,0] = numt2-(aa-1)*(k) + PP[k,1] = bb*(k) + + if kkmid: + DP[k-1,0] = PP[k,0] - PC[0] + DP[k-1,1] = PP[k,1] - PC[1] + + kmax = int((PP[kmid,0] + PP[kmid-1,0])/2 ) + kmin = int((PP[kmid,0] + PP[kmid+1,0])/2 ) + return [kmin,kmax,PP,DP] + +def SpreadPoint3D(M,R,sdep): + + [row,col,dep,numt2] = M.shape + PS = np.zeros([row,col,dep,numt2],dtype=complex) + inx = 0 + iny = 0 + + for k in range(np.mod(inx,R),row,R): + for ss in range(np.mod(iny,R),col,R): + PS[k,ss,:,:] = 1 + iny = iny + 1 + inx = inx + 1 + + + for k in range(numt2): + PS[:,:,:,k] = np.fft.ifftn(PS[:,:,:,k]) + PS[:,:,:,k] = np.fft.fftshift(PS[:,:,:,k]) + + + SPREAD = PS[:,:,sdep,:] + SPREAD = np.fft.ifft(SPREAD,axis=2) + SPREAD = np.fft.fftshift(SPREAD,axes=2) + + return SPREAD + +def SpreadPoint(M,R,scol): + + [row,col,numt2] = M.shape + PS = np.zeros([row,col,numt2],dtype=complex) + inx = 0 + + for l in range(0,numt2): + for k in range(np.mod(inx,R),row,R): + PS[k,:,l] = 1 + inx = inx + 1 + + #PS = 1*(M!=0) + #PS = 0*M + 1 + #SPREAD = KTT(PS,0) + + for k in range(numt2): + PS[:,:,k] = np.fft.ifft2(PS[:,:,k]) + PS[:,:,k] = np.fft.fftshift(PS[:,:,k]) + + + SPREAD = PS[:,scol,:] + SPREAD = np.fft.ifft(SPREAD,axis=1) + SPREAD = np.fft.fftshift(SPREAD,axes=1) + + return SPREAD + +def ConvolSP(M1,M2): + M2 = M2[::-1,::-1] + M3 = sc.signal.convolve2d(M1,M2, boundary='wrap', mode='same') + return M3 + +def KTBLASTMETHOD_4D_kxky(ITOT,R,mode): + + + ################################################################### + # Training Stage # + ################################################################### + [row,col,dep,numt2] = ITOT.shape + # INPUT PARAMETERS + iteshort = 1 + tetest = int(dep/2) + numt = int(numt2) + Dyy = int(row*0.1) + rmid = int(row/2) + cmid = int(col/2) + + TKdata = np.zeros(ITOT.shape,dtype=complex) + UKdata = np.zeros(ITOT.shape,dtype=complex) + Kdata = np.zeros(ITOT.shape,dtype=complex) + Kdata_NEW0 = np.zeros(ITOT.shape,dtype=complex) + Kdata_NEW = np.zeros(ITOT.shape,dtype=complex) + KTBLAST0 = np.zeros(ITOT.shape,dtype=complex) + KTBLAST = np.zeros(ITOT.shape,dtype=complex) + + + for k in range(numt2): + # THE FULL KSPACE + Kdata[:,:,:,k] = np.fft.fftn(ITOT[:,:,:,k]) + Kdata[:,:,:,k] = np.fft.fftshift(Kdata[:,:,:,k]) + # UNDERSAMPLING STEP AND FILLED WITH ZEROS THE REST + UKdata[:,:,:,k] = UNDER(Kdata[:,:,:,k],mode,R,k) + + # GENERATING THE TRAINING DATA WITH SUBSAMPLING the Center IN KX , KY + for k in range(numt): + TKdata[rmid-Dyy:rmid+Dyy+1,cmid-Dyy:cmid+Dyy+1,:,k] = Kdata[rmid-Dyy:rmid+Dyy+1,cmid-Dyy:cmid+Dyy+1,:,k] + + [kmin,kmax,PP,DP] = get_points4D(row,col,dep,numt2,R,mode) + + if iteshort==1: + print(PP) + print(DP) + print('range of k = ',kmin,kmax) + + SPREAD = SpreadPoint3D(UKdata,R,tetest) + ################################################################### + # RECONSTRUCTION # + ################################################################### + ZE1 = iteshort + (tetest-1)*(iteshort) + ZE2 = (tetest+1)*(iteshort) + (dep)*(1-iteshort) + + for zi in range(ZE1,ZE2): + + if rank==0: + print('4D KTBLAST: R = ' + str(R) + ' and z = ' + str(zi)+'/'+str(dep)) + + ### CONSTRUCT THE REFERENCE M_TRAINING + B2 = KTT3D(TKdata,zi) + B2 = FilteringHigh(B2,0.3) + M2 = 4*np.abs(B2)**2 + #M2 = GetSymmetric(M2) + ### INTERPOLATE IF NUMT0.001) + + if scantype=='0G': + PHASE0[:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,k])>0.001) + 10*varPHASE0 + PHASE1[:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,k])>0.001) + 10*varPHASE0 + np.pi*Sq[:,:,k]/VENC + + if scantype=='-G+G': + PHASE0[:,:,k] = gamma*B0*TE*np.ones([row,col]) + 10*varPHASE0 - np.pi*Sq[:,:,k]/VENC + PHASE1[:,:,k] = gamma*B0*TE*np.ones([row,col]) + 10*varPHASE0 + np.pi*Sq[:,:,k]/VENC + + RHO0[:,:,k] = modulus*np.cos(PHASE0[:,:,k]) + Drho + 1j*modulus*np.sin(PHASE0[:,:,k]) + 1j*Drho2 + RHO1[:,:,k] = modulus*np.cos(PHASE1[:,:,k]) + Drho + 1j*modulus*np.sin(PHASE1[:,:,k]) + 1j*Drho2 + + + if np.ndim(Sq)==4: + [row,col,dep,numt2] = Sq.shape + [X,Y,Z] = np.meshgrid(np.linspace(0,col,col),np.linspace(0,row,row),np.linspace(0,dep,dep)) + + for k in range(numt2): + + if noise: + Drho = np.random.normal(0,0.2,[row,col,dep]) + Drho2 = np.random.normal(0,0.2,[row,col,dep]) + else: + Drho = np.zeros([row,col,dep]) + Drho2 = np.zeros([row,col,dep]) + + varPHASE0 = np.random.randint(-10,11,size=(row,col,dep))*np.pi/180*(np.abs(Sq[:,:,:,k])<0.001) + modulus = 0.5 + 0.5*(np.abs(Sq[:,:,:,k])>0.001) + + if scantype=='0G': + PHASE0[:,:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,:,k])>0.001) + 10*varPHASE0 + PHASE1[:,:,:,k] = (gamma*B0*TE+0.01*X)*(np.abs(Sq[:,:,:,k])>0.001) + 10*varPHASE0 + np.pi*Sq[:,:,:,k]/VENC + + if scantype=='-G+G': + PHASE0[:,:,:,k] = gamma*B0*TE*np.ones([row,col,dep]) + varPHASE0 - np.pi*Sq[:,:,:,k]/VENC + PHASE1[:,:,:,k] = gamma*B0*TE*np.ones([row,col,dep]) + varPHASE0 + np.pi*Sq[:,:,:,k]/VENC + + RHO0[:,:,:,k] = modulus*np.cos(PHASE0[:,:,:,k]) + Drho + 1j*modulus*np.sin(PHASE0[:,:,:,k]) + 1j*Drho2 + RHO1[:,:,:,k] = modulus*np.cos(PHASE1[:,:,:,k]) + Drho + 1j*modulus*np.sin(PHASE1[:,:,:,k]) + 1j*Drho2 + + + + + return [RHO0,RHO1] + +def undersampling(Sqx,Sqy,Sqz,options,savepath): + + R = options['kt-BLAST']['R'] + mode = options['kt-BLAST']['mode'] + transpose = True + + for r in R: + + if rank==0: + print('Using Acceleration Factor R = ' + str(r)) + print('Component x of M0') + + [M0,M1] = GenerateMagnetization(Sqx,options['kt-BLAST']['VENC'],options['kt-BLAST']['noise'],scantype='0G') + if transpose: + M0 = M0.transpose((0,2,1,3)) + M1 = M1.transpose((0,2,1,3)) + + if mode=='ky': + M0_kt = KTBLASTMETHOD_4D_ky(M0,r,mode) + if mode=='kxky': + M0_kt = KTBLASTMETHOD_4D_kxky(M0,r,mode) + + if rank==0: + print('\n Component x of M1') + + if mode=='ky': + M1_kt = KTBLASTMETHOD_4D_ky(M1,r,mode) + if mode=='kxky': + M1_kt = KTBLASTMETHOD_4D_kxky(M1,r,mode) + + + Sqx_kt = phase_contrast(M1_kt,M0_kt,options['kt-BLAST']['VENC'],scantype='0G') + + del M0,M1 + del M0_kt, M1_kt + + [M0,M1] = GenerateMagnetization(Sqy,options['kt-BLAST']['VENC'],options['kt-BLAST']['noise'],scantype='0G') + if transpose: + M0 = M0.transpose((0,2,1,3)) + M1 = M1.transpose((0,2,1,3)) + if rank==0: + print('\n Component y of M0') + + if mode=='ky': + M0_kt = KTBLASTMETHOD_4D_ky(M0,r,mode) + if mode=='kxky': + M0_kt = KTBLASTMETHOD_4D_kxky(M0,r,mode) + + + if rank==0: + print('\n Component y of M1') + + if mode=='ky': + M1_kt = KTBLASTMETHOD_4D_ky(M1,r,mode) + if mode=='kxky': + M1_kt = KTBLASTMETHOD_4D_kxky(M1,r,mode) + + Sqy_kt = phase_contrast(M1_kt,M0_kt,options['kt-BLAST']['VENC'],scantype='0G') + + del M0,M1 + del M0_kt, M1_kt + + [M0,M1] = GenerateMagnetization(Sqz,options['kt-BLAST']['VENC'],options['kt-BLAST']['noise'],scantype='0G') + if transpose: + M0 = M0.transpose((0,2,1,3)) + M1 = M1.transpose((0,2,1,3)) + if rank==0: + print('\n Component z of M0') + + if mode=='ky': + M0_kt = KTBLASTMETHOD_4D_ky(M0,r,mode) + + if mode=='kxky': + M0_kt = KTBLASTMETHOD_4D_kxky(M0,r,mode) + + if rank==0: + print('\n Component z of M1') + + if mode=='ky': + M1_kt = KTBLASTMETHOD_4D_ky(M1,r,mode) + if mode=='kxky': + M1_kt = KTBLASTMETHOD_4D_kxky(M1,r,mode) + + if rank==0: + print(' ') + + Sqz_kt = phase_contrast(M1_kt,M0_kt,options['kt-BLAST']['VENC'],scantype='0G') + + + if transpose: + Sqx_kt = Sqx_kt.transpose((0,2,1,3)) + Sqy_kt = Sqy_kt.transpose((0,2,1,3)) + Sqz_kt = Sqz_kt.transpose((0,2,1,3)) + + + if options['kt-BLAST']['save']: + if rank==0: + print('saving the sequences in ' + savepath) + seqname = options['kt-BLAST']['name'] +'_R' + str(r) + '.npz' + print('sequence name: ' + seqname) + np.savez_compressed( savepath + seqname, x=Sqx_kt, y=Sqy_kt,z=Sqz_kt) + + del Sqx_kt,Sqy_kt,Sqz_kt + + + + diff --git a/codes/mesh_generator.py b/codes/mesh_generator.py new file mode 100755 index 0000000..5207cc8 --- /dev/null +++ b/codes/mesh_generator.py @@ -0,0 +1,71 @@ +from dolfin import * + + +mesh_in = '/home/yeye/Desktop/leomesh.xml' +mesh_out = '/home/yeye/Desktop/aorta.h5' +mesh = Mesh(mesh_in) + +hdf = HDF5File(mesh.mpi_comm(), mesh_out, 'w') +boundaries = MeshFunction('size_t', mesh,2) + +marked = 1 +testmesh = 0 + +hdf.write(mesh, '/mesh') + +if marked==1: + + class Inlet(SubDomain): + def inside(self, x, on_boundary): + return on_boundary and between(x[0],(0.1975,0.1989)) and between(x[2],(0.07,0.1)) + + class Outlet(SubDomain): + def inside(self, x, on_boundary): + return on_boundary and between(x[0],(0.1975,0.1989)) and between(x[2],(0,0.04)) + + class Walls(SubDomain): + def inside(self, x, on_boundary): + return on_boundary + + + outlet = Outlet() + inlet = Inlet() + walls = Walls() + + boundaries.set_all(0) + walls.mark(boundaries,1) + outlet.mark(boundaries,3) + inlet.mark(boundaries,2) + + + + +hdf.write(boundaries, '/boundaries') +hdf.close() + + + +if testmesh: + print('Testing Mesh...') + meshname = mesh_out + pathtoB = '/home/yeye/Desktop/boundaries.xdmf' + mesh = Mesh() + hdf = HDF5File(mesh.mpi_comm(), meshname , 'r') + hdf.read(mesh, '/mesh', False) + boundaries = MeshFunction('size_t', mesh , mesh.topology().dim() - 1) + hdf.read(boundaries, '/boundaries') + # To save the boundaries information + XDMFFile(pathtoB).write(boundaries) + print('Boundary info printed in ' + pathtoB) + + + + + + + + + + + +