Overhauled to simplify creation of materials with averaged dielectric properties.

gprMax/yee_cell_build.pyx

@@ -22,17 +22,96 @@ from gprMax.materials import Material
 from gprMax.yee_cell_setget_rigid cimport get_rigid_Ex, get_rigid_Ey, get_rigid_Ez, get_rigid_Hx, get_rigid_Hy, get_rigid_Hz
 
 
+cpdef void create_electric_average(int i, int j, int k, int numID1, int numID2, int numID3, int numID4, int componentID, G):
+    """This function creates a new material by averaging the dielectric properties of the surrounding cells.
+        
+    Args:
+        i, j, k (int): Cell coordinates.
+        numID1, numID2, numID3, numID4 (int): Numeric IDs for materials in surrounding cells.
+        componentID (int): Numeric ID for electric field component.
+        G (class): Grid class instance - holds essential parameters describing the model.
+    """
+
+    # Make an ID composed of the names of the four materials that will be averaged
+    requiredID = G.materials[numID1].ID + '+' + G.materials[numID2].ID + '+' + G.materials[numID3].ID + '+' + G.materials[numID4].ID
+    
+    # Check if this material already exists
+    tmp = requiredID.split('+')
+    material = [x for x in G.materials if
+             x.ID.count(tmp[0]) == requiredID.count(tmp[0]) and
+             x.ID.count(tmp[1]) == requiredID.count(tmp[1]) and
+             x.ID.count(tmp[2]) == requiredID.count(tmp[2]) and
+             x.ID.count(tmp[3]) == requiredID.count(tmp[3])]
+             
+    if material:
+        G.ID[componentID, i, j, k] = material[0].numID
+    else:
+        # Create new material
+        newNumID = len(G.materials)
+        m = Material(newNumID, requiredID, G)
+        # Create averaged constituents for material
+        m.er = np.mean((G.materials[numID1].er, G.materials[numID2].er, G.materials[numID3].er, G.materials[numID4].er), axis=0)
+        m.se = np.mean((G.materials[numID1].se, G.materials[numID2].se, G.materials[numID3].se, G.materials[numID4].se), axis=0)
+        m.mr = np.mean((G.materials[numID1].mr, G.materials[numID2].mr, G.materials[numID3].mr, G.materials[numID4].mr), axis=0)
+        m.sm = np.mean((G.materials[numID1].sm, G.materials[numID2].sm, G.materials[numID3].sm, G.materials[numID4].sm), axis=0)
+
+        # Append the new material object to the materials list
+        G.materials.append(m)
+        
+        G.ID[componentID, i, j, k] = newNumID
+
+
+cpdef void create_magnetic_average(int i, int j, int k, int numID1, int numID2, int componentID, G):
+    """This function creates a new material by averaging the dielectric properties of the surrounding cells.
+        
+    Args:
+        i, j, k (int): Cell coordinates.
+        numID1, numID2 (int): Numeric IDs for materials in surrounding cells.
+        componentID (int): Numeric ID for electric field component.
+        G (class): Grid class instance - holds essential parameters describing the model.
+    """
+
+    # Make an ID composed of the names of the two materials that will be averaged
+    requiredID = G.materials[numID1].ID + '+' + G.materials[numID2].ID
+    
+    # Check if this material already exists
+    tmp = requiredID.split('+')
+    material = [x for x in G.materials if
+                (x.ID.count(tmp[0]) == requiredID.count(tmp[0]) and
+                 x.ID.count(tmp[1]) == requiredID.count(tmp[1])) or
+                (x.ID.count(tmp[0]) % 2 == 0 and x.ID.count(tmp[1]) % 2 == 0)]
+    
+    if material:
+        G.ID[componentID, i, j, k] = material[0].numID
+    else:
+        # Create new material
+        newNumID = len(G.materials)
+        m = Material(newNumID, requiredID, G)
+        # Create averaged constituents for material
+        m.er = np.mean((G.materials[numID1].er, G.materials[numID2].er), axis=0)
+        m.se = np.mean((G.materials[numID1].se, G.materials[numID2].se), axis=0)
+        m.mr = np.mean((G.materials[numID1].mr, G.materials[numID2].mr), axis=0)
+        m.sm = np.mean((G.materials[numID1].sm, G.materials[numID2].sm), axis=0)
+        
+        # Append the new material object to the materials list
+        G.materials.append(m)
+        
+        G.ID[componentID, i, j, k] = newNumID
+
+
 cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidE, np.uint32_t[:, :, :, ::1] ID, G):
     """This function builds the electric field components in the ID array.
         
     Args:
         solid, rigid, ID (memoryviews): Access to solid, rigid and ID arrays
+        G (class): Grid class instance - holds essential parameters describing the model.
     """
     
     cdef Py_ssize_t i, j, k
     cdef int numID1, numID2, numID3, numID4
 
     # Ex component
+    componentID = 0
     for i in range(0, G.nx):
         for j in range(1, G.ny):
             for k in range(1, G.nz):
@@ -48,36 +127,13 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                     
                     # If all values are the same no need to average
                     if numID1 == numID2 and numID1 == numID3 and numID1 == numID4:
-                        ID[0, i, j, k] = numID1
+                        ID[componentID, i, j, k] = numID1
                     else:
                         # Averaging is required
-                        # Make an ID composed of the names of the four materials that will be averaged
-                        requiredID = G.materials[numID1].ID + '+' + G.materials[numID2].ID + '+' + G.materials[numID3].ID + '+' + G.materials[numID4].ID
-                        # Check if this material already exists
-                        tmp = requiredID.split('+')
-                        material = [x for x in G.materials if
-                                 x.ID.count(tmp[0]) == requiredID.count(tmp[0]) and
-                                 x.ID.count(tmp[1]) == requiredID.count(tmp[1]) and
-                                 x.ID.count(tmp[2]) == requiredID.count(tmp[2]) and
-                                 x.ID.count(tmp[3]) == requiredID.count(tmp[3])]
-                                 
-                        if material:
-                            ID[0, i, j, k] = material[0].numID
-                        else:
-                            # Create new material
-                            newNumID = len(G.materials)
-                            m = Material(newNumID, requiredID, G)
-                            # Create averaged constituents for material
-                            m.er = np.mean((G.materials[numID1].er, G.materials[numID2].er, G.materials[numID3].er, G.materials[numID4].er), axis=0)
-                            m.se = np.mean((G.materials[numID1].se, G.materials[numID2].se, G.materials[numID3].se, G.materials[numID4].se), axis=0)
-                            m.mr = np.mean((G.materials[numID1].mr, G.materials[numID2].mr, G.materials[numID3].mr, G.materials[numID4].mr), axis=0)
-                            m.sm = np.mean((G.materials[numID1].sm, G.materials[numID2].sm, G.materials[numID3].sm, G.materials[numID4].sm), axis=0)
-                            
-                            # Append the new material object to the materials list
-                            G.materials.append(m)
-                            ID[0, i, j, k] = newNumID
+                        create_electric_average(i, j, k, numID1, numID2, numID3, numID4, componentID, G)
 
     # Ey component
+    componentID = 1
     for i in range(1, G.nx):
         for j in range(0, G.ny):
             for k in range(1, G.nz):
@@ -90,39 +146,16 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                     numID2 = solid[i - 1, j, k]
                     numID3 = solid[i - 1, j, k - 1]
                     numID4 = solid[i, j, k - 1]
-
+                    
                     # If all values are the same no need to average
                     if numID1 == numID2 and numID1 == numID3 and numID1 == numID4:
-                        ID[1, i, j, k] = numID1
+                        ID[componentID, i, j, k] = numID1
                     else:
                         # Averaging is required
-                        # Make an ID composed of the names of the four materials that will be averaged
-                        requiredID = G.materials[numID1].ID + '+' + G.materials[numID2].ID + '+' + G.materials[numID3].ID + '+' + G.materials[numID4].ID
-                        # Check if this material already exists
-                        tmp = requiredID.split('+')
-                        material = [x for x in G.materials if
-                                 x.ID.count(tmp[0]) == requiredID.count(tmp[0]) and
-                                 x.ID.count(tmp[1]) == requiredID.count(tmp[1]) and
-                                 x.ID.count(tmp[2]) == requiredID.count(tmp[2]) and
-                                 x.ID.count(tmp[3]) == requiredID.count(tmp[3])]
-                                 
-                        if material:
-                            ID[1, i, j, k] = material[0].numID
-                        else:
-                            # Create new material
-                            newNumID = len(G.materials)
-                            m = Material(newNumID, requiredID, G)
-                            # Create averaged constituents for material
-                            m.er = np.mean((G.materials[numID1].er, G.materials[numID2].er, G.materials[numID3].er, G.materials[numID4].er), axis=0)
-                            m.se = np.mean((G.materials[numID1].se, G.materials[numID2].se, G.materials[numID3].se, G.materials[numID4].se), axis=0)
-                            m.mr = np.mean((G.materials[numID1].mr, G.materials[numID2].mr, G.materials[numID3].mr, G.materials[numID4].mr), axis=0)
-                            m.sm = np.mean((G.materials[numID1].sm, G.materials[numID2].sm, G.materials[numID3].sm, G.materials[numID4].sm), axis=0)
-                            
-                            # Append the new material object to the materials list
-                            G.materials.append(m)
-                            ID[1, i, j, k] = newNumID
+                        create_electric_average(i, j, k, numID1, numID2, numID3, numID4, componentID, G)
 
     # Ez component
+    componentID = 2
     for i in range(1, G.nx):
         for j in range(1, G.ny):
             for k in range(0, G.nz):
@@ -138,34 +171,10 @@ cpdef void build_electric_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                     
                     # If all values are the same no need to average
                     if numID1 == numID2 and numID1 == numID3 and numID1 == numID4:
-                        ID[2, i, j, k] = numID1
+                        ID[componentID, i, j, k] = numID1
                     else:
                         # Averaging is required
-                        # Make an ID composed of the names of the four materials that will be averaged
-                        requiredID = G.materials[numID1].ID + '+' + G.materials[numID2].ID + '+' + G.materials[numID3].ID + '+' + G.materials[numID4].ID
-                        # Check if this material already exists
-                        tmp = requiredID.split('+')
-                        material = [x for x in G.materials if
-                                 x.ID.count(tmp[0]) == requiredID.count(tmp[0]) and
-                                 x.ID.count(tmp[1]) == requiredID.count(tmp[1]) and
-                                 x.ID.count(tmp[2]) == requiredID.count(tmp[2]) and
-                                 x.ID.count(tmp[3]) == requiredID.count(tmp[3])]
-                                 
-                        if material:
-                            ID[2, i, j, k] = material[0].numID
-                        else:
-                            # Create new material
-                            newNumID = len(G.materials)
-                            m = Material(newNumID, requiredID, G)
-                            # Create averaged constituents for material
-                            m.er = np.mean((G.materials[numID1].er, G.materials[numID2].er, G.materials[numID3].er, G.materials[numID4].er), axis=0)
-                            m.se = np.mean((G.materials[numID1].se, G.materials[numID2].se, G.materials[numID3].se, G.materials[numID4].se), axis=0)
-                            m.mr = np.mean((G.materials[numID1].mr, G.materials[numID2].mr, G.materials[numID3].mr, G.materials[numID4].mr), axis=0)
-                            m.sm = np.mean((G.materials[numID1].sm, G.materials[numID2].sm, G.materials[numID3].sm, G.materials[numID4].sm), axis=0)
-
-                            # Append the new material object to the materials list
-                            G.materials.append(m)
-                            ID[2, i, j, k] = newNumID
+                        create_electric_average(i, j, k, numID1, numID2, numID3, numID4, componentID, G)
 
 
 cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:, :, :, ::1] rigidH, np.uint32_t[:, :, :, ::1] ID, G):
@@ -179,6 +188,7 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
     cdef int numID1, numID2
 
     # Hx component
+    componentID = 3
     for i in range(1, G.nx):
         for j in range(0, G.ny):
             for k in range(0, G.nz):
@@ -192,35 +202,13 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                     
                     # If all values are the same no need to average
                     if numID1 == numID2:
-                        ID[3, i, j, k] = numID1
+                        ID[componentID, i, j, k] = numID1
                     else:
                         # Averaging is required
-                        # Make an ID composed of the names of the two materials that will be averaged
-                        requiredID = G.materials[numID1].ID + '+' + G.materials[numID2].ID
-                        # Check if this material already exists
-                        tmp = requiredID.split('+')
-                        material = [x for x in G.materials if
-                                    (x.ID.count(tmp[0]) == requiredID.count(tmp[0]) and
-                                     x.ID.count(tmp[1]) == requiredID.count(tmp[1])) or
-                                    (x.ID.count(tmp[0]) % 2 == 0 and x.ID.count(tmp[1]) % 2 == 0)]
-                        
-                        if material:
-                            ID[3, i, j, k] = material[0].numID
-                        else:
-                            # Create new material
-                            newNumID = len(G.materials)
-                            m = Material(newNumID, requiredID, G)
-                            # Create averaged constituents for material
-                            m.er = np.mean((G.materials[numID1].er, G.materials[numID2].er), axis=0)
-                            m.se = np.mean((G.materials[numID1].se, G.materials[numID2].se), axis=0)
-                            m.mr = np.mean((G.materials[numID1].mr, G.materials[numID2].mr), axis=0)
-                            m.sm = np.mean((G.materials[numID1].sm, G.materials[numID2].sm), axis=0)
-                            
-                            # Append the new material object to the materials list
-                            G.materials.append(m)
-                            ID[3, i, j, k] = newNumID
+                        create_magnetic_average(i, j, k, numID1, numID2, componentID, G)
 
     # Hy component
+    componentID = 4
     for i in range(0, G.nx):
         for j in range(1, G.ny):
             for k in range(0, G.nz):
@@ -237,32 +225,10 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                         ID[4, i, j, k] = numID1
                     else:
                         # Averaging is required
-                        # Make an ID composed of the names of the two materials that will be averaged
-                        requiredID = G.materials[numID1].ID + '+' + G.materials[numID2].ID
-                        # Check if this material already exists
-                        tmp = requiredID.split('+')
-                        material = [x for x in G.materials if
-                                    (x.ID.count(tmp[0]) == requiredID.count(tmp[0]) and
-                                     x.ID.count(tmp[1]) == requiredID.count(tmp[1])) or
-                                    (x.ID.count(tmp[0]) % 2 == 0 and x.ID.count(tmp[1]) % 2 == 0)]
-                        
-                        if material:
-                            ID[4, i, j, k] = material[0].numID
-                        else:
-                            # Create new material
-                            newNumID = len(G.materials)
-                            m = Material(newNumID, requiredID, G)
-                            # Create averaged constituents for material
-                            m.er = np.mean((G.materials[numID1].er, G.materials[numID2].er), axis=0)
-                            m.se = np.mean((G.materials[numID1].se, G.materials[numID2].se), axis=0)
-                            m.mr = np.mean((G.materials[numID1].mr, G.materials[numID2].mr), axis=0)
-                            m.sm = np.mean((G.materials[numID1].sm, G.materials[numID2].sm), axis=0)
-                            
-                            # Append the new material object to the materials list
-                            G.materials.append(m)
-                            ID[4, i, j, k] = newNumID
+                        create_magnetic_average(i, j, k, numID1, numID2, componentID, G)
 
     # Hz component
+    componentID = 5
     for i in range(0, G.nx):
         for j in range(0, G.ny):
             for k in range(1, G.nz):
@@ -279,27 +245,5 @@ cpdef void build_magnetic_components(np.uint32_t[:, :, ::1] solid, np.int8_t[:,
                         ID[5, i, j, k] = numID1
                     else:
                         # Averaging is required
-                        # Make an ID composed of the names of the two materials that will be averaged
-                        requiredID = G.materials[numID1].ID + '+' + G.materials[numID2].ID
-                        # Check if this material already exists
-                        tmp = requiredID.split('+')
-                        material = [x for x in G.materials if
-                                    (x.ID.count(tmp[0]) == requiredID.count(tmp[0]) and
-                                     x.ID.count(tmp[1]) == requiredID.count(tmp[1])) or
-                                    (x.ID.count(tmp[0]) % 2 == 0 and x.ID.count(tmp[1]) % 2 == 0)]
-                        
-                        if material:
-                            ID[5, i, j, k] = material[0].numID
-                        else:
-                            # Create new material
-                            newNumID = len(G.materials)
-                            m = Material(newNumID, requiredID, G)
-                            # Create averaged constituents for material
-                            m.er = np.mean((G.materials[numID1].er, G.materials[numID2].er), axis=0)
-                            m.se = np.mean((G.materials[numID1].se, G.materials[numID2].se), axis=0)
-                            m.mr = np.mean((G.materials[numID1].mr, G.materials[numID2].mr), axis=0)
-                            m.sm = np.mean((G.materials[numID1].sm, G.materials[numID2].sm), axis=0)
-                            
-                            # Append the new material object to the materials list
-                            G.materials.append(m)
-                            ID[5, i, j, k] = newNumID
+                        create_magnetic_average(i, j, k, numID1, numID2, componentID, G)
+