Assymetric pattern fix

.gitignore

@@ -1,2 +1,4 @@
 access-token.txt
 local-build.sh
+OSS_platform/**/__pycache__
+

OSS_platform/Axon_files/axon.py

@@ -138,21 +138,21 @@ def get_axonparams(self):
             'total_nodes'    int
                 Total number of compartments
             'ranvier_nodes'  int
-                Number of node of Ranvier compartments
+                Number of node of Ranvier compartments 
             'para1_nodes'    int
-                Number of first paranodal compartments
+                Number of first paranodal compartments (MYSA)
             'para2_nodes'    int
-                Number of second paranodal compartments
+                Number of second paranodal compartments (FLUT)
             'inter_nodes'    int
-                Number of internodal compartments
+                Number of internodal compartments (STIN)
             'ranvier_length' float
                 Length of the Node of Ranvier
             'para1_length'   float
-                Length of the first paranodal compartments
+                Length of the first paranodal compartments (MYSA)
             'para2_length'   float
-                Length of the second paranodal compartments
+                Length of the second paranodal compartments (FLUT)
             'inter_length'   float
-                Length of the internodal compartments
+                Length of the internodal compartments (STIN)
             'deltax'         float
                 Length from a node of Ranvier to the next
             'fiberD'         float
@@ -162,9 +162,9 @@ def get_axonparams(self):
             'node_diameter': float
                 Diameter of the nodes of Ranvier
             'para1_diameter': float
-                Diameter of the first paranodal compartments
+                Diameter of the first paranodal compartments (MYSA)
             'para2_diameter': float,
-                Diameter of the second paranodal compartments
+                Diameter of the second paranodal compartments (FLUT)
             'condg':          float
                 Axial conductivity
             'lamellas':       int
@@ -191,23 +191,40 @@ def __create_ref_nodes(self):
         ncomp = nstin+nmysa+nflut+1      
 
         ref_nodes = np.zeros(n,)
-        for i in range(0, n):
-            j = i % ncomp
-            if j == 0 or j == 1:  # mysa <-> ranvier node
-                if i == 0:  # first of all nodes
-                    ref_nodes[i] = l_ranvier/2.
-                else:
-                    ref_nodes[i] = ref_nodes[i-1]+l_para1/2.+l_ranvier/2.
-            elif (j > 1 and j < (int(nmysa/2)+1)) or (j > (int(nmysa/2)+nflut+nstin) and j <(nmysa+nflut+nstin)): # mysa <-> mysa node
-                ref_nodes[i] = ref_nodes[i-1]+l_para1
-            elif j == (int(nmysa/2)+1) or j == (int(nmysa/2)+nflut+nstin):  # mysa <-> flut node
-                ref_nodes[i] = ref_nodes[i-1]+l_para1/2.+l_para2/2.
-            elif (j > (int(nmysa/2)+1) and j < (int(nmysa/2)+int(nflut/2)+1)) or (j > (int(nmysa/2)+int(nflut/2)+nstin) and j < (int(nmysa/2)+nflut+nstin)): # flut <-> flut node
-                ref_nodes[i] = ref_nodes[i-1]+l_para2
-            elif j == (int(nmysa/2)+int(nflut/2)+1) or j == (int(nmysa/2)+int(nflut/2)+nstin):  # flut <-> stin node
-                ref_nodes[i] = ref_nodes[i-1]+l_para2/2.+l_inter/2.
-            else:  # stin <-> stin node
-                ref_nodes[i] = ref_nodes[i-1]+l_inter
+
+        # An easier implementation based on the structure of the internodals segments
+        structure = 'RMF'+6*'S'+'FM' # Node of Rnavier + internodal structure
+
+        # Let's rename the lengths, just for convenience
+        L_R = l_ranvier 
+        L_F = l_para2
+        L_M = l_para1
+        L_S = l_inter
+
+        compartment = [0] # intial node of Ranvier 
+        for i in range(n-1):
+            l1 = eval('L_'+structure[i%len(structure)]) # current segment length
+            l2 = eval('L_'+structure[(i+1)%len(structure)]) # next segment length
+            compartment.append(compartment[-1]+(l1+l2)/2) 
+
+        ref_nodes = np.array(compartment)
+        # for i in range(0, n):
+        #     j = i % ncomp
+        #     if j == 0 or j == 1:  # mysa <-> ranvier node
+        #         if i == 0:  # first of all nodes
+        #             ref_nodes[i] = l_ranvier/2.
+        #         else:
+        #             ref_nodes[i] = ref_nodes[i-1]+l_para1/2.+l_ranvier/2.
+        #     elif (j > 1 and j < (int(nmysa/2)+1)) or (j > (int(nmysa/2)+nflut+nstin) and j <(nmysa+nflut+nstin)): # mysa <-> mysa node
+        #         ref_nodes[i] = ref_nodes[i-1]+l_para1
+        #     elif j == (int(nmysa/2)+1) or j == (int(nmysa/2)+nflut+nstin):  # mysa <-> flut node
+        #         ref_nodes[i] = ref_nodes[i-1]+l_para1/2.+l_para2/2.
+        #     elif (j > (int(nmysa/2)+1) and j < (int(nmysa/2)+int(nflut/2)+1)) or (j > (int(nmysa/2)+int(nflut/2)+nstin) and j < (int(nmysa/2)+nflut+nstin)): # flut <-> flut node
+        #         ref_nodes[i] = ref_nodes[i-1]+l_para2
+        #     elif j == (int(nmysa/2)+int(nflut/2)+1) or j == (int(nmysa/2)+int(nflut/2)+nstin):  # flut <-> stin node
+        #         ref_nodes[i] = ref_nodes[i-1]+l_para2/2.+l_inter/2.
+        #     else:  # stin <-> stin node
+        #         ref_nodes[i] = ref_nodes[i-1]+l_inter
 
         return ref_nodes * 1e-6  # um -> m
 
@@ -217,6 +234,9 @@ def __create_nodes(self):
                                                        self.__ref_nodes[-1]/2
         else:
             ref_nodes = self.__ref_nodes
+
+        # The reference nodes are not 100% symmetric, i.e. sometimes ref_nodes[i]!=ref_nodes[-i]
+        # due to the fitting option
         # standard position(along x axis in [x,y,z])
         nodes = np.zeros((len(ref_nodes),3))
         nodes[:,0] = ref_nodes
@@ -256,17 +276,17 @@ def __create_axon_parameters(self, diameter):
 
         axon_parameters = copy.deepcopy(self.AXONPARAMETERS[diameter])
 
-        nranvier = axon_parameters['ranvier_nodes']
-        nstin = int(float(axon_parameters['inter_nodes'])/(nranvier-1))
-        nmysa = int(float(axon_parameters['para1_nodes'])/(nranvier-1))
-        nflut = int(float(axon_parameters['para2_nodes'])/(nranvier-1))
+        nranvier = axon_parameters['ranvier_nodes'] #21
+        nstin = int(float(axon_parameters['inter_nodes'])/(nranvier-1)) # 6
+        nmysa = int(float(axon_parameters['para1_nodes'])/(nranvier-1)) # 2
+        nflut = int(float(axon_parameters['para2_nodes'])/(nranvier-1)) # 2
 
         axon_parameters['fiberD']=diameter
         axon_parameters['total_nodes'] = axon_parameters['ranvier_nodes'] + \
             axon_parameters['para1_nodes']+axon_parameters['para2_nodes'] + \
-            axon_parameters['inter_nodes']
+            axon_parameters['inter_nodes'] #221
         axon_parameters['inter_length'] = (axon_parameters['deltax'] -
                                        axon_parameters['ranvier_length'] -
                                        nmysa*axon_parameters['para1_length'] -
-                                       nflut*axon_parameters['para2_length'])/float(nstin)
-        return axon_parameters
+                                       nflut*axon_parameters['para2_length'])/float(nstin) #70.5
+        return axon_parameters

OSS_platform/CSF_refinement_new.py

@@ -304,6 +304,7 @@ def Refine_CSF(MRI_param,DTI_param,Scaling,Domains,Field_calc_param,rel_div,CSF_
         z_neuron_min=min(Vertices_neur[:,2])
 
         #go over all voxels and check whether it contains CSF and intersect with the mesh
+        eps = 0.000000001
         for x_coord in x_vect:
             for y_coord in y_vect:
                 for z_coord in z_vect:
@@ -314,9 +315,9 @@ def Refine_CSF(MRI_param,DTI_param,Scaling,Domains,Field_calc_param,rel_div,CSF_
 
                     if (x_pos<=x_neuron_max+CSF_ref_add and x_pos>=x_neuron_min-CSF_ref_add and y_pos<=y_neuron_max+CSF_ref_add and y_pos>=y_neuron_min-CSF_ref_add and z_pos<=z_neuron_max+CSF_ref_add and z_pos>=z_neuron_min-CSF_ref_add):
 
-                        xv_mri=int((x_coord)/MRI_param.x_vox_size-0.000000001)                                  #defines number of steps to get to the voxels containing x[0] coordinate
-                        yv_mri=(int((y_coord)/MRI_param.y_vox_size-0.000000001))*MRI_param.M_x                  #defines number of steps to get to the voxels containing x[0] and x[1] coordinates
-                        zv_mri=(int((z_coord)/MRI_param.z_vox_size-0.000000001))*MRI_param.M_x*MRI_param.M_y          #defines number of steps to get to the voxels containing x[0], x[1] and x[2] coordinates
+                        xv_mri= x_coord//(MRI_param.x_vox_size-eps)                                 #defines number of steps to get to the voxels containing x[0] coordinate
+                        yv_mri=(y_coord//(MRI_param.y_vox_size-eps))*MRI_param.M_x                  #defines number of steps to get to the voxels containing x[0] and x[1] coordinates
+                        zv_mri=(z_coord//(MRI_param.z_vox_size-eps))*MRI_param.M_x*MRI_param.M_y    #defines number of steps to get to the voxels containing x[0], x[1] and x[2] coordinates
                         glob_index=int(xv_mri+yv_mri+zv_mri)
 
                         pnt=Point(x_pos,y_pos,z_pos)

OSS_platform/Dict_corrector.py

@@ -4,17 +4,17 @@
 
 @author: konstantin
 """
+import numpy as np
 
 #converts from user friendly units to the platform's format
-
+    
 def rearrange_Inp_dict(dict_from_GUI):
 
     #fitst, change empty lines to 0
     for key in dict_from_GUI:
         if dict_from_GUI[key]=='' or dict_from_GUI[key]=='0':
             dict_from_GUI[key]=0 
 
-    import numpy as np
     #second, some angles will be changed to rads
     if dict_from_GUI['Global_rot'] == 1:
         dict_from_GUI['alpha_array_glob'] =[i*np.pi/(180.0) for i in dict_from_GUI['alpha_array_glob']]

OSS_platform/FEM_in_spectrum.py

@@ -198,8 +198,12 @@ class Conductivity : public dolfin::Expression
 
     if int(sine_freq)==int(signal_freq):        
         c_unscaled = CompiledExpression(compile_cpp_code(conductivity_code).Conductivity(),
-                       c00=c_unscaled00, c01=c_unscaled01, c02=c_unscaled02, c11=c_unscaled11, c12=c_unscaled12, c22=c_unscaled22, degree=0)
-        tensor= as_tensor([[c_unscaled[0], c_unscaled[1], c_unscaled[2]], [c_unscaled[1], c_unscaled[3], c_unscaled[4]],[c_unscaled[2],c_unscaled[4],c_unscaled[5]]])
+                                        c00=c_unscaled00, c01=c_unscaled01, c02=c_unscaled02,
+                                        c11=c_unscaled11, c12=c_unscaled12, c22=c_unscaled22,
+                                        degree=0)
+        tensor = as_tensor([[c_unscaled[0], c_unscaled[1], c_unscaled[2]], 
+                            [c_unscaled[1], c_unscaled[3], c_unscaled[4]],
+                            [c_unscaled[2],c_unscaled[4],c_unscaled[5]]])
         f_vector_repr=project(tensor,TensorFunctionSpace(mesh, "Lagrange", 1),solver_type="cg", preconditioner_type="amg")
         file=File('Tensors/Ellipsoids_unscaled_at_'+str(signal_freq)+'_Hz.pvd')
         file<<f_vector_repr    
@@ -423,8 +427,16 @@ def solve_Laplace(Sim_setup,Solver_type,Vertices_array,Domains,core,Full_IFFT,ou
     Sim_setup.mesh.coordinates()
     Sim_setup.mesh.init()
 
-    # to get conductivity (and permittivity if EQS formulation) mapped accrodingly to the subdomains. k_val_r is just a list of conductivities (S/mm!) in a specific order to scale the cond. tensor
-    kappa,k_val_r=get_dielectric_properties_from_subdomains(Sim_setup.mesh,Sim_setup.subdomains,Sim_setup.Laplace_eq,Domains.Float_contacts,Sim_setup.conductivities,Sim_setup.rel_permittivities,Sim_setup.sine_freq)
+    # to get conductivity (and permittivity if EQS formulation) mapped accrodingly to the subdomains. 
+    # k_val_r is just a list of conductivities (S/mm!) in a specific order to scale the cond. tensor
+    kappa,k_val_r = get_dielectric_properties_from_subdomains(Sim_setup.mesh, 
+                                                              Sim_setup.subdomains,
+                                                              Sim_setup.Laplace_eq,
+                                                              Domains.Float_contacts,
+                                                              Sim_setup.conductivities,
+                                                              Sim_setup.rel_permittivities,
+                                                              Sim_setup.sine_freq)
+
     if int(Sim_setup.sine_freq)==int(Sim_setup.signal_freq):
         file=File('Field_solutions/Conductivity_map_'+str(Sim_setup.signal_freq)+'Hz.pvd')
         file<<kappa[0]
@@ -434,19 +446,33 @@ def solve_Laplace(Sim_setup,Solver_type,Vertices_array,Domains,core,Full_IFFT,ou
 
     # to get tensor scaled by the conductivity map
     if Sim_setup.anisotropy==1:
-        Cond_tensor=get_scaled_cond_tensor(Sim_setup.mesh,Sim_setup.subdomains,Sim_setup.sine_freq,Sim_setup.signal_freq,Sim_setup.unscaled_tensor,k_val_r)
+        Cond_tensor = get_scaled_cond_tensor(Sim_setup.mesh,
+                                             Sim_setup.subdomains,
+                                             Sim_setup.sine_freq,
+                                             Sim_setup.signal_freq,
+                                             Sim_setup.unscaled_tensor,
+                                             k_val_r)
     else:
-        Cond_tensor=False  #just to initialize
-
-    #In case of current-controlled stimulation, Dirichlet_bc or the whole potential distribution will be scaled afterwards (due to the system's linearity)
-    V_space,Dirichlet_bc,ground_index=get_solution_space_and_Dirichlet_BC(Sim_setup.c_c,Sim_setup.mesh,Sim_setup.boundaries,Sim_setup.element_order,Sim_setup.Laplace_eq,Domains.Contacts,Domains.fi)
+        Cond_tensor = False  #just to initialize
+
+    # In case of current-controlled stimulation, Dirichlet_bc 
+    # the whole potential distribution will be scaled afterwards 
+    # (due to the system's linearity)
+    V_space, Dirichlet_bc, ground_index = get_solution_space_and_Dirichlet_BC(Sim_setup.c_c,
+                                                                              Sim_setup.mesh,
+                                                                              Sim_setup.boundaries,
+                                                                              Sim_setup.element_order,
+                                                                              Sim_setup.Laplace_eq,
+                                                                              Domains.Contacts,
+                                                                              Domains.fi)
     #ground index refers to the ground in .med/.msh file
 
     # to solve the Laplace equation div(kappa*grad(phi))=0   (variational form: a(u,v)=L(v))    
-    phi_sol=define_variational_form_and_solve(V_space,Dirichlet_bc,kappa,Sim_setup.Laplace_eq,Cond_tensor,Solver_type)        
+    phi_sol =define_variational_form_and_solve(V_space, Dirichlet_bc, kappa, Sim_setup.Laplace_eq, 
+                                               Cond_tensor, Solver_type)        
 
     if Sim_setup.Laplace_eq=='EQS':
-        (phi_r,phi_i)=phi_sol.split(deepcopy=True)
+        (phi_r,phi_i)= phi_sol.split(deepcopy=True)
     else:
         phi_r=phi_sol
         phi_i=Function(V_space)
@@ -460,7 +486,8 @@ def solve_Laplace(Sim_setup,Solver_type,Vertices_array,Domains,core,Full_IFFT,ou
 
 
     if Sim_setup.c_c==1 or Sim_setup.CPE_status==1:     #we compute E-field, currents and impedances only for current-controlled or if CPE is used
-        J_ground=get_current(Sim_setup.mesh,Sim_setup.boundaries,Sim_setup.element_order,Sim_setup.Laplace_eq,Domains.Contacts,kappa,Cond_tensor,phi_r,phi_i,ground_index)                
+        J_ground = get_current(Sim_setup.mesh, Sim_setup.boundaries, Sim_setup.element_order, Sim_setup.Laplace_eq,
+                               Domains.Contacts, kappa, Cond_tensor, phi_r, phi_i, ground_index)            
         #If EQS, J_ground is a complex number
 
         #V_across=max(Domains.fi[:], key=abs)        # voltage drop in the system
@@ -484,25 +511,34 @@ def solve_Laplace(Sim_setup,Solver_type,Vertices_array,Domains,core,Full_IFFT,ou
                 print("Currently, CPE can be used only for simulations with two contacts. Please, assign the rest to 'None'")
                 raise SystemExit
 
-            Dirichlet_bc_with_CPE,total_impedance=get_CPE_corrected_Dirichlet_BC(Sim_setup.boundaries,Sim_setup.CPE_param,Sim_setup.Laplace_eq,Sim_setup.sine_freq,Sim_setup.signal_freq,Domains.Contacts,Domains.fi,V_across,Z_tissue,V_space)
+            Dirichlet_bc_with_CPE, total_impedance = get_CPE_corrected_Dirichlet_BC(Sim_setup.boundaries,
+                                                                                    Sim_setup.CPE_param,
+                                                                                    Sim_setup.Laplace_eq,
+                                                                                    Sim_setup.sine_freq,
+                                                                                    Sim_setup.signal_freq,
+                                                                                    Domains.Contacts, Domains.fi,
+                                                                                    V_across, Z_tissue, V_space)
 
             f=open('Field_solutions/Impedance'+str(core)+'.csv','ab')
             np.savetxt(f, total_impedance, delimiter=" ")
             f.close()
 
             # to solve the Laplace equation for the adjusted Dirichlet   
-            phi_sol_CPE=define_variational_form_and_solve(V_space,Dirichlet_bc_with_CPE,kappa,Sim_setup.Laplace_eq,Cond_tensor,Solver_type)   
+            phi_sol_CPE = define_variational_form_and_solve(V_space, Dirichlet_bc_with_CPE, kappa,
+                                                            Sim_setup.Laplace_eq, Cond_tensor, Solver_type)   
             if Sim_setup.Laplace_eq=='EQS':
                 (phi_r_CPE,phi_i_CPE)=phi_sol_CPE.split(deepcopy=True)
             else:
                 phi_r_CPE=phi_sol_CPE
                 phi_i_CPE=Function(V_space)
                 phi_i_CPE.vector()[:] = 0.0
 
-            J_ground_CPE=get_current(Sim_setup.mesh,Sim_setup.boundaries,Sim_setup.element_order,Sim_setup.Laplace_eq,Domains.Contacts,kappa,Cond_tensor,phi_r_CPE,phi_i_CPE,ground_index)                
+            J_ground_CPE = get_current(Sim_setup.mesh, Sim_setup.boundaries, Sim_setup.element_order, 
+                                       Sim_setup.Laplace_eq, Domains.Contacts, kappa, Cond_tensor, 
+                                       phi_r_CPE, phi_i_CPE, ground_index)                
 
             # just resaving
-            phi_sol,phi_r,phi_i,J_ground=(phi_sol_CPE,phi_r_CPE,phi_i_CPE,J_ground_CPE)
+            phi_sol, phi_r, phi_i, J_ground = (phi_sol_CPE, phi_r_CPE, phi_i_CPE, J_ground_CPE)
 
     if Full_IFFT==1:
         Hdf=HDF5File(Sim_setup.mesh.mpi_comm(), "Field_solutions_functions/solution"+str(np.round(Sim_setup.sine_freq,6))+".h5", "w")
@@ -567,7 +603,8 @@ def solve_Laplace(Sim_setup,Solver_type,Vertices_array,Domains,core,Full_IFFT,ou
                 else:
                     Dirichlet_bc_scaled.append(DirichletBC(V_space, np.real((Domains.fi[bc_i])/J_ground),Sim_setup.boundaries,Domains.Contacts[bc_i]))
 
-        phi_sol_check=define_variational_form_and_solve(V_space,Dirichlet_bc_scaled,kappa,Sim_setup.Laplace_eq,Cond_tensor,Solver_type)        
+        phi_sol_check = define_variational_form_and_solve(V_space, Dirichlet_bc_scaled, kappa, 
+                                                          Sim_setup.Laplace_eq, Cond_tensor, Solver_type)        
 
         if Sim_setup.Laplace_eq=='EQS':
             (phi_r_check,phi_i_check)=phi_sol_check.split(deepcopy=True)