From 1346cbc8259a033649a7cdac44eba4098c60fd53 Mon Sep 17 00:00:00 2001 From: Craig Warren Date: Thu, 25 Feb 2016 18:38:57 +0000 Subject: [PATCH] Overhauled to simplify creation of materials with averaged dielectric properties. --- gprMax/yee_cell_build.pyx | 248 +++++++++++++++----------------------- 1 file changed, 96 insertions(+), 152 deletions(-) diff --git a/gprMax/yee_cell_build.pyx b/gprMax/yee_cell_build.pyx index 82e2dfcd..9739df10 100644 --- a/gprMax/yee_cell_build.pyx +++ b/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) +