diff --git a/lbmpy/sparse/mapping.py b/lbmpy/sparse/mapping.py index 3cff6fd21e98d1c91df48b848a8b821fd4dd85dd..a43aba664a6efddd813577c24664cb7b1c095a45 100644 --- a/lbmpy/sparse/mapping.py +++ b/lbmpy/sparse/mapping.py @@ -108,11 +108,7 @@ class SparseLbMapper: # Add fluid cells coordinates_fluid = np.argwhere(np.bitwise_and(self._flag_arr, self.fluid_flag)).astype(np.uint32) - #print("coordinates_fluid", coordinates_fluid) - #print("self._flag_arr", self._flag_arr, " boundary_mask", self.boundary_mask, " density:", self.density_flag, " ubb_fag", self.ubb_flag) - #print("bitw_and _flag_arr, boundary_mask:", np.bitwise_and(self._flag_arr, self.boundary_mask)) coordinates_boundary = np.argwhere(np.bitwise_and(self._flag_arr, self.boundary_mask)).astype(np.uint32) - #print("Coordinates boundary:", coordinates_boundary) self._num_fluid_cells = coordinates_fluid.shape[0] total_cells = len(coordinates_fluid) + len(coordinates_boundary) @@ -131,21 +127,14 @@ class SparseLbMapper: flag_arr = self.flag_array no_slip_flag = self.no_slip_flag fluid_boundary_mask = self.ubb_flag | self.fluid_flag | self.density_flag - #print("Fuid boundary mask", fluid_boundary_mask) - #fluid_boundary_mask = self.fluid_flag - #print("Flag_arr.shape", flag_arr.shape) result = [] print("flag:", self.flag_array) for direction_idx, direction in enumerate(stencil): - #print("Direction:", direction_idx, direction) if all(d_i == 0 for d_i in direction): assert direction_idx == 0 continue for own_cell_idx, cell in enumerate(self.fluid_coordinates): domain_size = (len(self.flag_array), len(self.flag_array[0])) - #if (direction[0] != 0 and (domain_size[0]-1)*(1-direction[0])/2 == cell[0]) or (direction[1] != 0 and (domain_size[1]-1)*(1-direction[1])/2 == cell[1]): - # print("direction:", direction_idx, "cell: ", cell) - # continue test = [] inv_neighbor_cell = np.array([cell_i - dir_i for cell_i, dir_i in zip(cell, direction)]) test.append(("Zelle", own_cell_idx)) @@ -161,15 +150,11 @@ class SparseLbMapper: at_border = False for i, x_i in enumerate(inv_neighbor_cell): if x_i == (ghost_layers - 1): #if inv_neighbor_cell NOT noslip/fuid & one coordinate 0 - #print("before", inv_neighbor_cell[i]) inv_neighbor_cell[i] += flag_arr.shape[i] - (2 * ghost_layers) - #print("after", inv_neighbor_cell[i]) at_border = True test.append(("At Border Else", neighbor_cell_idx, "First")) elif x_i == flag_arr.shape[i] - ghost_layers: - #print("before", inv_neighbor_cell[i]) inv_neighbor_cell[i] -= flag_arr.shape[i] - (2 * ghost_layers) - #print("after", inv_neighbor_cell[i]) at_border = True test.append(("At Border Else", neighbor_cell_idx, "Second")) if at_border: @@ -203,88 +188,78 @@ class SparseLbPeriodicityMapper: raise IndexError("Coordinate not found") else: return self.mapping._sorter[left] - + def create_index_arr(self): stencil = self.method.stencil print("domain_size:", self.domain_size) result = [] + inner = [] for direction_idx, direction in enumerate(stencil): - if all(d_i == 0 for d_i in direction) or direction_idx != 5: - #assert direction_idx == 0 + if all(d_i == 0 for d_i in direction): continue - print("direction:", direction, ", ", direction_idx) - for pos in range(0, 2): + #print("direction:", direction, ", ", direction_idx) + for pos in range(0,2): sop = (pos+1)%2 if direction[pos] != 0: + # periodic/parallel index_array = [] - print("break") coord = int((self.domain_size[pos]-1)*(1-direction[pos])/2) prev_read = [0,0] for i in range(0, self.domain_size[sop]): - write = [1,0] + write = [0,0] cur_read = [0,0] write[pos] = coord write[sop] = i - cur_read[pos] = (write[pos]-direction[pos])%self.domain_size[pos] - cur_read[sop] = (write[sop]-direction[sop])%self.domain_size[sop] + cur_read = [(write_i - dir_i)%ds_i for write_i, dir_i, ds_i in zip(write, direction, self.domain_size)] if cur_read[pos] < prev_read[pos] or cur_read[sop] < prev_read[sop]: result.append(tuple(index_array)) index_array = [] - print("break") - print("write:", write, "read:", cur_read) + #print("write:", write, "read:", cur_read) write_idx = self.cell_idx(tuple(write)) read_idx = self.cell_idx(tuple(cur_read)) index_array.append((direction_idx, pdf_index(write_idx, direction_idx, len(self.mapping)), pdf_index(read_idx, direction_idx, len(self.mapping)))) + #index_array.append([direction_idx, write, cur_read]) prev_read[pos] = cur_read[pos] prev_read[sop] = cur_read[sop] result.append(tuple(index_array)) + # inner + pos_bound = int((self.domain_size[pos]-1)*(1+direction[pos])/2) + pos_mid = pos_bound+direction[pos]*(-self.domain_size[pos]+1) + sop_position = [int((self.domain_size[sop]-1)*(1+direction[sop])/2)] if direction[sop] != 0 else [0, self.domain_size[sop]-1] + #print("pos_bound:", pos_bound, "pos_mid", pos_mid, "sop_position:", sop_position) + for b in sop_position: + for i in range(pos_bound, pos_mid, -direction[pos]): + write = [0,0] + write[pos] = i + write[sop] = b + read = [write_i - dir_i for write_i, dir_i in zip(write, direction)] + #print("write:", write, "read:", read) + inner.append([direction_idx, write, read]) + if direction[pos] == 0: + # inner + pos_low = 1 + pos_high = self.domain_size[pos]-1 + sop_position = int((self.domain_size[sop]-1)*(1+direction[sop])/2) + #print("pos_low:", pos_low, "pos_high:", pos_high, "sop_position:", sop_position) + for i in range(pos_low, pos_high): + write = [0,0] + write[pos] = i + write[sop] = sop_position + read = [write_i - dir_i for write_i, dir_i in zip(write, direction)] + #print("write:", write, "read:", read) + inner.append([direction_idx, write, read]) result = list(dict.fromkeys(result)) list_result = [] - print("as tuples") - for i_a in result: - print(i_a) + for index_array in result: list_index_array = [] - for ass in i_a: - list_index_array.append(list(ass)) + for write_read_pair in index_array: + list_index_array.append(list(write_read_pair)) list_result.append(list_index_array) - print("as lists") + list_result.append(inner) for i_a in list_result: print(i_a) return list_result - - - """ - def create_index_arr(self): - stencil = self.method.stencil - print("domain_size:", self.domain_size) - result = [] - for direction_idx, direction in enumerate(stencil): - if all(d_i == 0 for d_i in direction): - assert direction_idx == 0 - continue - print("direction:", direction, ", ", direction_idx) - for own_cell_idx, cell in enumerate(self.mapping.fluid_coordinates): - from_cell_idx = ((cell[0]-direction[0])%self.domain_size[0], (cell[1]-direction[1])%self.domain_size[1]) - write = pdf_index(own_cell_idx, direction_idx, len(self.mapping)) - read = pdf_index(self.cell_idx(from_cell_idx), direction_idx, len(self.mapping)) - if direction[1] == 1 and cell[1] == 0: - print("(own_cell_idx:", own_cell_idx, "), to cell:", cell, "from cell:", from_cell_idx) - print("write:", write, "read:", read) - result.append([direction_idx, write, read]) - elif direction[1] == -1 and cell[1] == self.domain_size[1]-1: - print("(own_cell_idx:", own_cell_idx, "), to cell:", cell, "from cell:", from_cell_idx) - print("write:", write, "read:", read) - result.append([direction_idx, write, read]) - elif direction[0] == 1 and cell[0] == 0: - print("(own_cell_idx:", own_cell_idx, "), to cell:", cell, "from cell:", from_cell_idx) - print("write:", write, "read:", read) - result.append([direction_idx, write, read]) - elif direction[0] == -1 and cell[0] == self.domain_size[0]-1: - print("(own_cell_idx:", own_cell_idx, "), to cell:", cell, "from cell:", from_cell_idx) - print("write:", write, "read:", read) - result.append([direction_idx, write, read]) - """ - + class SparseLbBoundaryMapper: NEIGHBOR_IDX_NAME = 'nidx{}' diff --git a/lbmpy_tests/test_sparse_lbm.ipynb b/lbmpy_tests/test_sparse_lbm.ipynb index a12e549a1779d583c7a1c8fb8ee311ff73e9299b..cbbbbda9bb1542917565fd5677d376560f5abeed 100644 --- a/lbmpy_tests/test_sparse_lbm.ipynb +++ b/lbmpy_tests/test_sparse_lbm.ipynb @@ -53,7 +53,7 @@ "metadata": {}, "outputs": [], "source": [ - "domain_size = (2,3)\n", + "domain_size = (4,3)\n", "omega = 1.8\n", "target = 'cpu'\n", "\n", @@ -142,7 +142,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAFpCAYAAACBJomJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dbaxlV3kn+P+Tiq2eACObtjEe22ArKgUcOhim5DjiC2mSjG1FMURBsj+A5XZUIcLdINHSuIkUMtNfUCYQQcO45AwlbIkY0QMOJaYSx7EYOZF4ceEUfndTcQi+uNov0G2DSMdTVc98uLuYM5dbdap89r63bp3fT9o6Z++19l7r+oiXx89az67uDgAAAIv7qc2eAAAAwOlCgAUAADASARYAAMBIBFgAAAAjEWABAACMRIAFAAAwkoUDrKq6qKq+XFWPVtXDVfW+dfq8taqer6r9w/H7i44LAABwVFXtrqpnquqhY7S/rqq+UlX/VFX/dk3blVX1eFUdqKqbZ66/sqrurqpvDZ9nz5vHGBmsQ0k+0N2vT3JFkvdW1aXr9Pvr7r5sOP7XEcYFAAA46tNJrjxO+/eT/JskfzR7saq2JflkkquSXJrkupl45uYk93T39iT3DOfHtXCA1d0Hu/v+4fsPkjya5IJFnwsAAHCiuvverAZRx2p/prvvS/L/rGm6PMmB7n6iu19M8tkk1wxt1yS5bfh+W5K3z5vHqHuwquriJG9K8rV1mn+pqr5ZVX9eVT8/5rgAAAAv0QVJnpw5X8n/lzA6r7sPJquJpSSvmvewnx5rVlX18iSfT/L+7n5hTfP9SV7b3T+sqquT/FmS7cd4zs4kO5PkZS972f/4ute9bqwpAqeQB595erOnAKe9f/Gq8zZ7CsBEvvGNbzzX3edu9jxO1P/0yy/r733/8ELP+MYD//Rwkv82c+nW7r51oYeuqnWu9Ut92CgBVlWdkdXg6jPd/YW17bMBV3fvrar/varO6e7n1ul7a5Jbk2THjh29b9++MaYInGIu+dhHNnsKcNrb974PbPYUgIlU1T9s9hxOxnPfP5yv3XXhQs844/y/+2/dvWOkKc1aSXLRzPmFSZ4avj9dVed398GqOj/JM/MeNkYVwUryqSSPdvdHj9Hn1UO/VNXlw7jfW3RsAACABd2XZHtVXVJVZya5NsmeoW1PkuuH79cn+eK8h42RwXpLknclebCq9g/XPpjkNUnS3buS/FaS362qQ0n+Mcm13f2S024AAMBW0jncRyYdoaruSPLWJOdU1UqSDyU5I1mNSarq1Un2JfnvkxypqvcnubS7X6iqm5LclWRbkt3d/fDw2A8n+VxV3ZjkO0neOW8eCwdY3f03WX/d4myfTyT5xKJjAQAAW08nOfLStzWd2Bjd181p/89ZXf63XtveJHvXuf69JG87mXmMVuQCAADgWI5k2gzWqWLUMu0AAADLTAYLAACYVKdzeElKMAiwAACAyU29B+tUIcACAAAm1UkOC7AAAADGsSwZLEUuAAAARiKDBQAATKoTRS4AAADGshxvwRJgAQAAE+u0IhcAAACj6OTwcsRXilwAAACMRQYLAACYVMceLAAAgJFUDqc2exIbQoAFAABMqpMcsQcLAACAkyGDBQAATM4SQQAAgBF0BFgAAACjOdICLAAAgIUtUwZLkQsAAICRyGABAACT6lQOL0luR4AFAABMzh4sAACAESzTHiwBFgAAMLHK4V6OJYLL8VcCAABsABksAABgUp3kyJLkdgRYAADA5OzBAgAAGEG3PVgAAACcJBksAABgckcsEQQAAFjc6nuwlmPx3HL8lQAAwCZa3YO1yDF3hKrdVfVMVT10jPaqqo9X1YGqeqCq3jxc/7mq2j9zvFBV7x/a/qCqvjvTdvW8echgAQAAk9qgMu2fTvKJJLcfo/2qJNuH4xeT3JLkF7v78SSXJUlVbUvy3SR3ztz3x939Ryc6CRksAABgy+vue5N8/zhdrklye6/6apKzqur8NX3eluTvuvsfXuo8BFgAAMDkDnctdCQ5p6r2zRw7T3IKFyR5cuZ8Zbg269okd6y5dtOwpHB3VZ09bxBLBAEAgEl1aowiF891944F7l+vjGH/uLHqzCS/keTfzbTfkuTfD/3+fZKPJPlXxxtEgAUAAEzuyOa/aHglyUUz5xcmeWrm/Kok93f300cvzH6vqj9J8qV5g2z6XwkAAJzejpZpX+QYwZ4k7x6qCV6R5PnuPjjTfl3WLA9cs0frHUnWrVA4SwYLAADY8qrqjiRvzeperZUkH0pyRpJ0964ke5NcneRAkh8luWHm3p9J8qtJfmfNY/+wqi7Laoz47XXaf4IACwAAmFTnx4Uqphuj+7o57Z3kvcdo+1GSf77O9Xed7DwEWAAAwOQ24D1YpwQBFgAAMKnu5PDmF7nYEMvxVwIAAGwAGSwAAGBilSPrvobq9CPAAgAAJtWxRPCEVdVFVfXlqnq0qh6uqvet06eq6uNVdaCqHqiqNy86LgAAsHWcAu/B2hBjZLAOJflAd99fVa9I8o2quru7H5npc1WS7cPxi0luGT4BAIDTXKdyZOIy7aeKhUPB7j7Y3fcP33+Q5NEkF6zpdk2S23vVV5OcteatyAAAAFveqHuwquriJG9K8rU1TRckeXLmfGW4dnCdZ+xMsjNJXvOa14w5PQAAYJNspWV+ixjtr6yqlyf5fJL3d/cLa5vXuaXXe05339rdO7p7x7nnnjvW9AAAgE3SSY70Ty10bBWjZLCq6oysBlef6e4vrNNlJclFM+cXJnlqjLEBAIBTXeXwkpRpH6OKYCX5VJJHu/ujx+i2J8m7h2qCVyR5vrt/YnkgAABw+pHBOjlvSfKuJA9W1f7h2geTvCZJuntXkr1Jrk5yIMmPktwwwrgAAACnlIUDrO7+m6y/x2q2Tyd576JjAQAAW9OyLBEctYogAADAWt21pZb5LUKABQAATO7wkgRYy/FXAgAAbAAZLAAAYFKd5Ig9WAAAAGOopVkiKMACAAAmtfoeLBksAACAURxekvIPy/FXAgAAbAAZLAAAYFKdskQQAABgLEeWZPGcAAsAAJhUd3JYBgsAAGAcy7JEcDnydAAAABtABgsAAJjUapGL5cjtCLAAAIDJHc5yLBEUYAEAAJPq2IMFAADASZLBAgAAJrY8e7CW468EAAA21ZHUQsc8VbW7qp6pqoeO0V5V9fGqOlBVD1TVm2favl1VD1bV/qraN3P9lVV1d1V9a/g8e948BFgAAMCkjr5oeJHjBHw6yZXHab8qyfbh2JnkljXtv9zdl3X3jplrNye5p7u3J7lnOD8uARYAADC5I/1TCx3zdPe9Sb5/nC7XJLm9V301yVlVdf6cx16T5Lbh+21J3j5vHgIsAABgKzinqvbNHDtP8v4Lkjw5c74yXEtWCx3+ZVV9Y81zz+vug0kyfL5q3iCKXAAAAJNafdHwwmXan1uzfO9krTeBHj7f0t1PVdWrktxdVY8NGbGTJoMFAABMbuoiFydgJclFM+cXJnkqSbr76OczSe5McvnQ5+mjywiHz2fmDSLAAgAAJnX0RcOLHCPYk+TdQzXBK5I8390Hq+plVfWKJKmqlyX5tSQPzdxz/fD9+iRfnDeIJYIAAMCWV1V3JHlrVvdqrST5UJIzkqS7dyXZm+TqJAeS/CjJDcOt5yW5s6qS1fjoT7v7L4a2Dyf5XFXdmOQ7Sd45bx4CLAAAYHJTv2i4u6+b095J3rvO9SeSvPEY93wvydtOZh4CLAAAYFrjLfM75QmwAACASXUyVqGKU54ACwAAmNyyZLBUEQQAABiJDBYAADCpo2Xal4EACwAAmJwACwAAYAQdVQQBAABGsyxVBBW5AAAAGIkMFgAAMK22BwsAAGAUqggCAACMaFkCLHuwAAAARiKDBQAATEqZdgAAgBG1AAsAAGAcy/IeLAEWAAAwqV6iMu2KXAAAAIxklACrqnZX1TNV9dAx2t9aVc9X1f7h+P0xxgUAALaG7lro2CrGWiL46SSfSHL7cfr8dXf/+kjjAQAAW4Yqgielu++tqovHeBYAAHD62UpZqEVs5B6sX6qqb1bVn1fVzx+rU1XtrKp9VbXv2Wef3cDpAQAAU+isFrlY5NgqNirAuj/Ja7v7jUn+Q5I/O1bH7r61u3d0945zzz13g6YHAACwuA0JsLr7he7+4fB9b5IzquqcjRgbAADYZL1aqn2RY6vYkPdgVdWrkzzd3V1Vl2c1sPveRowNAABsPi8aPglVdUeStyY5p6pWknwoyRlJ0t27kvxWkt+tqkNJ/jHJtd1bKQ4FAABeqs7yFLkYq4rgdXPaP5HVMu4AAACnrQ1ZIggAACyzrVUJcBECLAAAYHLLskFIgAUAAEzOHiwAAIARrJZaX44Aa6NeNAwAAHDaE2ABAACTO9K10DFPVe2uqmeq6qFjtFdVfbyqDlTVA1X15uH6RVX15ap6tKoerqr3zdzzB1X13araPxxXz5uHAAsAAJjc6jLBl36cgE8nufI47Vcl2T4cO5PcMlw/lOQD3f36JFckeW9VXTpz3x9392XDsXfeJOzBAgAAJjf1HqzuvreqLj5Ol2uS3N7dneSrVXVWVZ3f3QeTHBye8YOqejTJBUkeeSnzkMECAAAm1al0L3YkOaeq9s0cO09yGhckeXLmfGW49mNDgPamJF+buXzTsKRwd1WdPW8QARYAALAVPNfdO2aOW0/y/vVSaD9efFhVL0/y+STv7+4Xhsu3JPnZJJdlNcv1kXmDWCIIAABM7hR4z/BKkotmzi9M8lSSVNUZWQ2uPtPdXzjaobufPvq9qv4kyZfmDSKDBQAATGt4D9aCSwQXtSfJu4dqglckeb67D1ZVJflUkke7+6OzN1TV+TOn70iyboXCWTJYAADA9CZOYVXVHUnemtW9WitJPpTkjCTp7l1J9ia5OsmBJD9KcsNw61uSvCvJg1W1f7j2waFi4B9W1WXD7L+d5HfmzUOABQAAbHndfd2c9k7y3nWu/03W35+V7n7Xyc5DgAUAAExu6jLtpwoBFgAAMLkTfFnwlifAAgAAJtWRwQIAABhHJ1mSAEuZdgAAgJHIYAEAAJOzBwsAAGAsAiwAAIAxlCIXAAAAo1mSDJYiFwAAACORwQIAAKbV3oMFAAAwniVZIijAAgAANsByZLDswQIAABiJDBYAADA9SwQBAABGIsACAAAYQSdRRRAAAGAcvSQZLEUuAAAARiKDBQAATG9JMlgCLAAAYHr2YAEAAIyjZLAAAABG0FmaJYKKXAAAAIxEBgsAAJhY2YMFAAAwmiVZIijAAgAAprckAZY9WAAAACORwQIAAKa3JBksARYAADCtjiIXAAAAY1mWFw2PsgerqnZX1TNV9dAx2quqPl5VB6rqgap68xjjAgAAW0QveMyxSExSVVdW1eND280z119ZVXdX1beGz7PnzWOsIhefTnLlcdqvSrJ9OHYmuWWkcQEAAJKXGJNU1bYknxzaL01yXVVdOtxzc5J7unt7knuG8+MaJcDq7nuTfP84Xa5Jcnuv+mqSs6rq/DHGBgAAWCAmuTzJge5+ortfTPLZoe/Re24bvt+W5O3z5rFRZdovSPLkzPnKcA0AAFgC1YsdIzhWTHK8WOW87j6YJMPnq+YNslFFLtYrGbLuP6aq2pnVlF1e85rXTDknAABgoyxeRfCcqto3c35rd996EvcfKyY54VjlRGxUgLWS5KKZ8wuTPLVex+Ef0q1JsmPHjiWpNQIAAKexEyxUMcdz3b1jgfuPFZOceYzrSfJ0VZ3f3QeH5YTPzBtko5YI7kny7qFyxxVJnj+aagMAANgAx4pJ7kuyvaouqaozk1w79D16z/XD9+uTfHHeIKNksKrqjiRvzWrabiXJh5KckSTdvSvJ3iRXJzmQ5EdJbhhjXAAAYIuYeG3aS41JuvtQVd2U5K4k25Ls7u6Hh8d+OMnnqurGJN9J8s558xglwOru6+a0d5L3jjEWAACw9Uz9ouFFYpLu3pvVAGzt9e8ledvJzGOj9mABAADLbEmqK2zUHiwAAIDTngwWAAAwvSXJYAmwAACASY34suBTngALAACY3uIvGt4SBFgAAMD0liSDpcgFAADASGSwAACAydmDBQAAMBYBFgAAwAiWqIqgPVgAAAAjkcECAACmtyQZLAEWAAAwPQEWAADAOOzBAgAA4KQIsAAAAEZiiSAAADC9JVkiKMACAACmtUTvwRJgAQAA0xNgAQAAjGRJAixFLgAAAEYigwUAAEyqYg8WAADAeARYAAAAI1iiKoL2YAEAAIxEBgsAAJjekmSwBFgAAMD0BFgAAADjWJY9WAIsAABgeksSYClyAQAAMBIZLAAAYFqdpclgCbAAAIDJLcseLEsEAQCA6fWCxxxVdWVVPV5VB6rq5nXaz66qO6vqgar6elW9Ybj+c1W1f+Z4oareP7T9QVV9d6bt6nnzkMECAAAmN2UGq6q2Jflkkl9NspLkvqra092PzHT7YJL93f2Oqnrd0P9t3f14kstmnvPdJHfO3PfH3f1HJzoXGSwAAGCruzzJge5+ortfTPLZJNes6XNpknuSpLsfS3JxVZ23ps/bkvxdd//DS52IAAsAAJje4ksEz6mqfTPHzpmnX5DkyZnzleHarG8m+c0kqarLk7w2yYVr+lyb5I41124alhXurqqz5/2ZAiwAAGBaiwZXqwHWc929Y+a4dWaEOsaosz6c5Oyq2p/kXyf52ySHfvyAqjOT/EaS/zhzzy1JfjarSwgPJvnIvD/VHiwAAGBSlfUjoBGtJLlo5vzCJE/NdujuF5LckCRVVUn+fjiOuirJ/d399Mw9P/5eVX+S5EvzJiKDBQAAbHX3JdleVZcMmahrk+yZ7VBVZw1tSfLbSe4dgq6jrsua5YFVdf7M6TuSPDRvIjJYAADA9CasItjdh6rqpiR3JdmWZHd3P1xV7xnadyV5fZLbq+pwkkeS3Hj0/qr6maxWIPydNY/+w6q6bJj9t9dp/wkCLAAAYHJTv2i4u/cm2bvm2q6Z719Jsv0Y9/4oyT9f5/q7TnYeAiwAAGB6EwdYpwoBFgAAML0lCbAUuQAAABiJDBYAADCtnn4P1qlCgAUAAExvSQKsUZYIVtWVVfV4VR2oqpvXaX9rVT1fVfuH4/fHGBcAANgaqhc7toqFM1hVtS3JJ7NaN34lyX1Vtae7H1nT9a+7+9cXHQ8AANiCtlCQtIgxMliXJznQ3U9094tJPpvkmhGeCwAAsKWMEWBdkOTJmfOV4dpav1RV36yqP6+qnz/Ww6pqZ1Xtq6p9zz777AjTAwAANtuyLBEcI8Cqda6t/Udwf5LXdvcbk/yHJH92rId1963dvaO7d5x77rkjTA8AANhUPcKxRYwRYK0kuWjm/MIkT8126O4XuvuHw/e9Sc6oqnNGGBsAANgKBFgn7L4k26vqkqo6M8m1SfbMdqiqV1dVDd8vH8b93ghjAwAAnDIWriLY3Yeq6qYkdyXZlmR3dz9cVe8Z2ncl+a0kv1tVh5L8Y5Jru3sLxaEAAMBLVdla+6gWMcqLhodlf3vXXNs18/0TST4xxlgAAMAWJMACAAAYRy3JAjYBFgAAMK0tVqhiEWMUuQAAACAyWAAAwAZQ5AIAAGAsAiwAAIBxyGABAACMZUkCLEUuAAAARiKDBQAATKstEQQAABiPAAsAAGBxleXJYNmDBQAAMBIZLAAAYHq9HCksARYAADC5ZVkiKMACAACm1VHkAgAAYCx1ZLNnsDEUuQAAABiJDBYAADC9JVkiKIMFAABMrnqxY+7zq66sqser6kBV3bxO+9lVdWdVPVBVX6+qN8y0fbuqHqyq/VW1b+b6K6vq7qr61vB59rx5CLAAAIBpdVbLtC9yHEdVbUvyySRXJbk0yXVVdemabh9Msr+7fyHJu5N8bE37L3f3Zd29Y+bazUnu6e7tSe4Zzo9LgAUAAExu4gzW5UkOdPcT3f1iks8muWZNn0uzGiSlux9LcnFVnTfnudckuW34fluSt8+biAALAADYCs6pqn0zx86ZtguSPDlzvjJcm/XNJL+ZJFV1eZLXJrlwaOskf1lV31jz3PO6+2CSDJ+vmjdJRS4AAIDpLV7k4rk1y/dm1QmM+OEkH6uq/UkeTPK3SQ4NbW/p7qeq6lVJ7q6qx7r73pcySQEWAAAwqcqJFapYwEqSi2bOL0zy1GyH7n4hyQ1JUlWV5O+HI9391PD5TFXdmdUlh/cmebqqzu/ug1V1fpJn5k3EEkEAAGBaixa4mFPkIsl9SbZX1SVVdWaSa5Psme1QVWcNbUny20nu7e4XquplVfWKoc/LkvxakoeGfnuSXD98vz7JF+dNRAYLAADY0rr7UFXdlOSuJNuS7O7uh6vqPUP7riSvT3J7VR1O8kiSG4fbz0ty52pSKz+d5E+7+y+Gtg8n+VxV3ZjkO0neOW8uAiwAAGByEy8RTHfvTbJ3zbVdM9+/kmT7Ovc9keSNx3jm95K87WTmIcACAACmN3GAdaoQYAEAAJObOoN1qhBgAQAA0+okR5YjwlJFEAAAYCQyWAAAwPSWI4ElwAIAAKZnDxYAAMBY5r8s+LQgwAIAACa3LBksRS4AAABGIoMFAABMq6PIBQAAwBgqSdmDBQAAMJIjmz2BjWEPFgAAwEhksAAAgMlZIggAADAGRS4AAADG0l40DAAAMBYvGgYAAOCkjBJgVdWVVfV4VR2oqpvXaa+q+vjQ/kBVvXmMcQEAgC2ie7Fji1h4iWBVbUvyySS/mmQlyX1Vtae7H5npdlWS7cPxi0luGT4BAIDTXSflPVgn7PIkB7r7ie5+Mclnk1yzps81SW7vVV9NclZVnT/C2AAAwFawJBmsMQKsC5I8OXO+Mlw72T5JkqraWVX7qmrfs88+O8L0AAAANsYYAVatc21tiHkifVYvdt/a3Tu6e8e555678OQAAIBTQC94bBFjlGlfSXLRzPmFSZ56CX0AAIDTVG2hZX6LGCODdV+S7VV1SVWdmeTaJHvW9NmT5N1DNcErkjzf3QdHGBsAANgKlmQP1sIZrO4+VFU3JbkrybYku7v74ap6z9C+K8neJFcnOZDkR0luWHRcAABgi+gkS1JFcIwlgunuvVkNomav7Zr53kneO8ZYAAAAp6pRAiwAAIBjqfTS7MESYAEAANMTYAEAAIxEgAUAADCCJSpyMUaZdgAAACLAAgAANkB1L3TMfX7VlVX1eFUdqKqb12k/u6rurKoHqurrVfWG4fpFVfXlqnq0qh6uqvfN3PMHVfXdqto/HFfPm4clggAAwPQm3INVVduSfDLJryZZSXJfVe3p7kdmun0wyf7ufkdVvW7o/7Ykh5J8oLvvr6pXJPlGVd09c+8fd/cfnehcZLAAAICJ9WqAtchxfJcnOdDdT3T3i0k+m+SaNX0uTXJPknT3Y0kurqrzuvtgd98/XP9BkkeTXPBS/1IBFgAAsBWcU1X7Zo6dM20XJHly5nwlPxkkfTPJbyZJVV2e5LVJLpztUFUXJ3lTkq/NXL5pWFa4u6rOnjdJARYAADCtzhgZrOe6e8fMcevMCHWMUWd9OMnZVbU/yb9O8rdZXR64+oCqlyf5fJL3d/cLw+VbkvxsksuSHEzykXl/qj1YAADA9KYt076S5KKZ8wuTPDXbYQiabkiSqqokfz8cqaozshpcfaa7vzBzz9NHv1fVnyT50ryJyGABAACTm7iK4H1JtlfVJVV1ZpJrk+z5/41fddbQliS/neTe7n5hCLY+leTR7v7omnvOnzl9R5KH5k1EBgsAAJjehFUEu/tQVd2U5K4k25Ls7u6Hq+o9Q/uuJK9PcntVHU7ySJIbh9vfkuRdSR4clg8myQe7e2+SP6yqy7K63PDbSX5n3lwEWAAAwJY3BER711zbNfP9K0m2r3Pf32T9PVzp7ned7DwEWAAAwLQ6yZHpMlinEgEWAAAwsRN6l9VpQYAFAABMT4AFAAAwkiUJsJRpBwAAGIkMFgAAMC1FLgAAAMbSSR/Z7ElsCAEWAAAwPXuwAAAAOBkyWAAAwLTswQIAABjRkiwRFGABAADTE2ABAACMoZcmwFLkAgAAYCQyWAAAwLQ6yRHvwQIAABjHkiwRFGABAADTE2ABAACMoZfmPViKXAAAAIxEBgsAAJhWJ92KXAAAAIxjSZYICrAAAIDpLUmRC3uwAAAARiKDBQAATKvbi4YBAABGsyRLBAVYAADA5FoGCwAAYAy9NBksRS4AAABGIoMFAABMq+M9WAAAAKNpe7AAAAAW1kl6STJYC+3BqqpXVtXdVfWt4fPsY/T7dlU9WFX7q2rfImMCAABbTPdqBmuRY46qurKqHq+qA1V18zrtZ1fVnVX1QFV9vareMO/eE413Zi1a5OLmJPd09/Yk9wznx/LL3X1Zd+9YcEwAAIAfq6ptST6Z5Koklya5rqouXdPtg0n2d/cvJHl3ko+dwL0nE+8kWTzAuibJbcP325K8fcHnAQAAp6E+0gsdc1ye5EB3P9HdLyb5bFZjlVmXZjVISnc/luTiqjpvzr0nHe8sGmCd190Hh0keTPKqY/TrJH9ZVd+oqp0LjgkAAGw10y4RvCDJkzPnK8O1Wd9M8ptJUlWXJ3ltkgvn3Hui8c6PzS1yUVV/leTV6zT93rx7Z7ylu5+qqlclubuqHuvue48x3s4kR4OwH1bV4ycxDvOdk+S5zZ4EJ8VvtrX4vbae0/Y3q/f/282ewlRO29/sNOX3msZrN3sCJ+MH+S93/VX/n+cs+Jh/tqaew63dfevwvdbpvzbt9eEkH6uq/UkeTPK3SQ6d4L0nbG6A1d2/cqy2qnq6qs7v7oNVdX6SZ47xjKeGz2eq6s6spuHWDbCGf0i3rtfG4qpqn31wW4vfbGvxe209frOtx2+2tfi9SJLuvnLiIVaSXDRzfmGSp9bM4YUkNyRJVVWSvx+OnznOvScU78xadIngniTXD9+vT/LFtR2q6mVV9Yqj35P8WpKHFhwXAADgqPuSbK+qS6rqzCTXZjVW+bGqOmtoS5LfTnLvEHQd79658c5ai74H68NJPldVNyb5TpJ3DpP/H5L8H919dZLzkty5GiTmp5P8aXf/xYLjAgAAJEm6+6ClGNsAAAMlSURBVFBV3ZTkriTbkuzu7oer6j1D+64kr09ye1UdTvJIkhuPd+/w6HXjneOp7uV44RerqmrnzFpVtgC/2dbi99p6/GZbj99sa/F7sWwEWAAAACNZdA8WAAAAAwHWEqqq/62qHquqB6rqzqo6a7PnxLFV1Tur6uGqOlJVqjCdwqrqyqp6vKoOVNXcN72zuapqd1U9U1UKL20BVXVRVX25qh4d/jvxfZs9J46vqv5ZVX29qr45/Gb/y2bPCTaCAGs53Z3kDd39C0n+U5J/t8nz4fgeyupL8dZ9tQGnhqraluSTSa7K6pvir6uqSzd3Vszx6SRTlw1mPIeSfKC7X5/kiiTv9Z+xU94/JfmX3f3GJJclubKqrtjkOcHkBFhLqLv/srsPDadfzWqtf05R3f1od3vh9qnv8iQHuvuJ7n4xyWeTXLPJc+I4hhfef3+z58GJ6e6D3X3/8P0HSR5NcsHmzorj6VU/HE7PGA6b/zntCbD4V0n+fLMnAaeBC5I8OXO+Ev/nDyZRVRcneVOSr23uTJinqrZV1f6svpz17u72m3HaW/Q9WJyiquqvkrx6nabf6+4vDn1+L6tLLj6zkXPjJ53I78Upr9a55t/Uwsiq6uVJPp/k/cMLQjmFdffhJJcN+73vrKo3dLd9j5zWBFinqe7+leO1V9X1SX49ydtarf5NN+/3YktYSXLRzPmFSZ7apLnAaamqzshqcPWZ7v7CZs+HE9fd/7Wq/u+s7nsUYHFas0RwCVXVlUn+5yS/0d0/2uz5wGniviTbq+qSqjozybVJ9mzynOC0UVWV5FNJHu3uj272fJivqs49Wqm4qv67JL+S5LHNnRVMT4C1nD6R5BVJ7q6q/VW1a7MnxLFV1TuqaiXJLyX5v6rqrs2eEz9pKBxzU5K7srr5/nPd/fDmzorjqao7knwlyc9V1UpV3bjZc+K43pLkXUn+5fC/Xfur6urNnhTHdX6SL1fVA1n9l1B3d/eXNnlOMLmyOgwAAGAcMlgAAAAjEWABAACMRIAFAAAwEgEWAADASARYAAAAIxFgAQAAjESABQAAMBIBFgAAwEj+X5LeU5OKFE0LAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAFpCAYAAACBJomJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3db6xlV3kn6N+biq1MgJZN2xiPy2ArKgUcpjF0yXHEF9IkGdsTxRAFyZYGLI+jChHuBolWy81IITP9BUWBDDSMLacpYasJiBlwKDGVOI6HkROJPy6bwn+xqDgEX1xtY5jBINJhquqdD3cXc/rm3jpVPnvfW7fO80hL5+y99t5rXR8BfnnXend1dwAAAFjcT231BAAAAM4UAiwAAICRCLAAAABGIsACAAAYiQALAABgJAIsAACAkSwcYFXVxVX1hap6vKoerap3rXPNG6vq+1V1cGi/t+i4AAAAx1XV3qp6tqoe2aD/VVX1xar6h6r612v6rqqqJ6rqUFXdMnP+pVV1T1V9Y/g8d948xshgHUnynu5+dZIrk7yzqi5b57q/6u7Lh/Y/jzAuAADAcR9PctUJ+r+X5F8l+cPZk1W1I8lHk1yd5LIk18/EM7ckube7dyW5dzg+oYUDrO4+3N0PDt9/kOTxJBct+lwAAICT1d33ZTWI2qj/2e6+P8n/u6briiSHuvvJ7v5xkk8luXbouzbJHcP3O5K8ed48Rt2DVVWXJHldki+v0/1LVfW1qvqzqvqFMccFAAB4gS5K8tTM8Ur+/4TRBd19OFlNLCV52byH/fRYs6qqFyf5TJJ3d/fza7ofTPLK7v5hVV2T5E+T7NrgOXuS7EmSF73oRf/8Va961VhT5Az28LPPbPUUAGAh/83LLtjqKbCNPPDAA8919/lbPY+T9d/+8ov6u987utAzHnjoHx5N8p9nTt3e3bcv9NBVtc65fqEPGyXAqqqzshpcfaK7P7u2fzbg6u79VfW/VtV53f3cOtfenuT2JNm9e3cfOHBgjClyhrv0Qx/Y6ikAwEIOvOs9Wz0FtpGq+rutnsOpeO57R/Plu3cu9IyzLvyb/9zdu0ea0qyVJBfPHO9M8vTw/ZmqurC7D1fVhUmenfewMaoIVpKPJXm8uz+4wTUvH65LVV0xjPvdRccGAABY0P1JdlXVpVV1dpLrkuwb+vYluWH4fkOSz8172BgZrDckeVuSh6vq4HDuvUlekSTdfVuS30ryu1V1JMnfJ7muu19w2g0AANhOOkf72KQjVNUnk7wxyXlVtZLkfUnOSlZjkqp6eZIDSf5JkmNV9e4kl3X381V1c5K7k+xIsre7Hx0e+/4kn66qm5J8K8lb581j4QCru/86669bnL3mI0k+suhYAADA9tNJjr3wbU0nN0b39XP6/1NWl/+t17c/yf51zn83yZtOZR6jFbkAAADYyLFMm8E6XYxaph0AAGCZyWABAACT6nSOLkkJBgEWAAAwuan3YJ0uBFgAAMCkOslRARYAAMA4liWDpcgFAADASGSwAACASXWiyAUAAMBYluMtWAIsAABgYp1W5AIAAGAUnRxdjvhKkQsAAICxyGABAACT6tiDBQAAMJLK0dRWT2JTCLAAAIBJdZJj9mABAABwKmSwAACAyVkiCAAAMIKOAAsAAGA0x1qABQAAsLBlymApcgEAADASGSwAAGBSncrRJcntCLAAAIDJ2YMFAAAwgmXagyXAAgAAJlY52suxRHA5/koAAIBNIIMFAABMqpMcW5LcjgALAACYnD1YAAAAI+i2BwsAAIBTJIMFAABM7pglggAAAItbfQ/WciyeW46/EgAA2EKre7AWaXNHqNpbVc9W1SMb9FdVfbiqDlXVQ1X1+uH8z1fVwZn2fFW9e+j7/ar69kzfNfPmIYMFAABMapPKtH88yUeS3LlB/9VJdg3tF5PcmuQXu/uJJJcnSVXtSPLtJHfN3PdH3f2HJzsJGSwAAGDb6+77knzvBJdcm+TOXvWlJOdU1YVrrnlTkr/p7r97ofMQYAEAAJM72rVQS3JeVR2YaXtOcQoXJXlq5nhlODfruiSfXHPu5mFJ4d6qOnfeIJYIAgAAk+rUGEUunuvu3Qvcv14Zw/5JZ9XZSX4jyb+d6b81yb8brvt3ST6Q5H840SACLAAAYHLHtv5FwytJLp453pnk6Znjq5M82N3PHD8x+72q/jjJ5+cNsuV/JQAAcGY7XqZ9kTaCfUnePlQTvDLJ97v78Ez/9VmzPHDNHq23JFm3QuEsGSwAAGDbq6pPJnljVvdqrSR5X5KzkqS7b0uyP8k1SQ4l+VGSG2fu/dkkv5rkd9Y89g+q6vKsxojfXKf/HxFgAQAAk+r8pFDFdGN0Xz+nv5O8c4O+HyX5p+ucf9upzkOABQAATG4T3oN1WhBgAQAAk+pOjm59kYtNsRx/JQAAwCaQwQIAACZWObbua6jOPAIsAABgUh1LBE9aVV1cVV+oqser6tGqetc611RVfbiqDlXVQ1X1+kXHBQAAto/T4D1Ym2KMDNaRJO/p7ger6iVJHqiqe7r7sZlrrk6ya2i/mOTW4RMAADjDdSrHJi7TfrpYOBTs7sPd/eDw/QdJHk9y0ZrLrk1yZ6/6UpJz1rwVGQAAYNsbdQ9WVV2S5HVJvrym66IkT80crwznDq/zjD1J9iTJK17xijGnBwAAbJHttMxvEaP9lVX14iSfSfLu7n5+bfc6t/R6z+nu27t7d3fvPv/888eaHgAAsEU6ybH+qYXadjFKBquqzspqcPWJ7v7sOpesJLl45nhnkqfHGBsAADjdVY4uSZn2MaoIVpKPJXm8uz+4wWX7krx9qCZ4ZZLvd/c/Wh4IAACceWSwTs0bkrwtycNVdXA4994kr0iS7r4tyf4k1yQ5lORHSW4cYVwAAIDTysIBVnf/ddbfYzV7TSd556JjAQAA29OyLBEctYogAADAWt21rZb5LUKABQAATO7okgRYy/FXAgAAbAIZLAAAYFKd5Jg9WAAAAGOopVkiKMACAAAmtfoeLBksAACAURxdkvIPy/FXAgAAbAIZLAAAYFKdskQQAABgLMeWZPGcAAsAAJhUd3JUBgsAAGAcy7JEcDnydAAAAJtABgsAAJjUapGL5cjtCLAAAIDJHc1yLBEUYAEAAJPq2IMFAADAKZLBAgAAJrY8e7CW468EAAC21LHUQm2eqtpbVc9W1SMb9FdVfbiqDlXVQ1X1+pm+b1bVw1V1sKoOzJx/aVXdU1XfGD7PnTcPARYAADCp4y8aXqSdhI8nueoE/Vcn2TW0PUluXdP/y919eXfvnjl3S5J7u3tXknuH4xMSYAEAAJM71j+1UJunu+9L8r0TXHJtkjt71ZeSnFNVF8557LVJ7hi+35HkzfPmIcACAAC2g/Oq6sBM23OK91+U5KmZ45XhXLJa6PAvquqBNc+9oLsPJ8nw+bJ5gyhyAQAATGr1RcMLl2l/bs3yvVO13gR6+HxDdz9dVS9Lck9VfX3IiJ0yGSwAAGByUxe5OAkrSS6eOd6Z5Okk6e7jn88muSvJFcM1zxxfRjh8PjtvEAEWAAAwqeMvGl6kjWBfkrcP1QSvTPL97j5cVS+qqpckSVW9KMmvJXlk5p4bhu83JPncvEEsEQQAALa9qvpkkjdmda/WSpL3JTkrSbr7tiT7k1yT5FCSHyW5cbj1giR3VVWyGh/9SXf/+dD3/iSfrqqbknwryVvnzUOABQAATG7qFw139/Vz+jvJO9c5/2SS125wz3eTvOlU5iHAAgAApjXeMr/TngALAACYVCdjFao47QmwAACAyS1LBksVQQAAgJHIYAEAAJM6XqZ9GQiwAACAyQmwAAAARtBRRRAAAGA0y1JFUJELAACAkchgAQAA02p7sAAAAEahiiAAAMCIliXAsgcLAABgJDJYAADApJRpBwAAGFELsAAAAMaxLO/BEmABAACT6iUq067IBQAAwEhGCbCqam9VPVtVj2zQ/8aq+n5VHRza740xLgAAsD1010JtuxhrieDHk3wkyZ0nuOavuvvXRxoPAADYNlQRPCXdfV9VXTLGswAAgDPPdspCLWIz92D9UlV9rar+rKp+YaOLqmpPVR2oqgPf+c53NnF6AADAFDqrRS4WadvFZgVYDyZ5ZXe/Nsm/T/KnG13Y3bd39+7u3n3++edv0vQAAAAWtykBVnc/390/HL7vT3JWVZ23GWMDAABbrFdLtS/StotNeQ9WVb08yTPd3VV1RVYDu+9uxtgAAMDW86LhU1BVn0zyxiTnVdVKkvclOStJuvu2JL+V5Her6kiSv09yXfd2ikMBAIAXqrM8RS7GqiJ4/Zz+j2S1jDsAAMAZa1OWCAIAAMtse1UCXIQACwAAmNyybBASYAEAAJOzBwsAAGAEq6XWlyPA2qwXDQMAAJzxBFgAAMDkjnUt1Oapqr1V9WxVPbJBf1XVh6vqUFU9VFWvH85fXFVfqKrHq+rRqnrXzD2/X1XfrqqDQ7tm3jwEWAAAwORWlwm+8HYSPp7kqhP0X51k19D2JLl1OH8kyXu6+9VJrkzyzqq6bOa+P+ruy4e2f94k7MECAAAmN/UerO6+r6ouOcEl1ya5s7s7yZeq6pyqurC7Dyc5PDzjB1X1eJKLkjz2QuYhgwUAAEyqU+lerCU5r6oOzLQ9pziNi5I8NXO8Mpz7iSFAe12SL8+cvnlYUri3qs6dN4gACwAA2A6e6+7dM+32U7x/vRTaTxYfVtWLk3wmybu7+/nh9K1Jfi7J5VnNcn1g3iCWCAIAAJM7Dd4zvJLk4pnjnUmeTpKqOiurwdUnuvuzxy/o7meOf6+qP07y+XmDyGABAADTGt6DteASwUXtS/L2oZrglUm+392Hq6qSfCzJ4939wdkbqurCmcO3JFm3QuEsGSwAAGB6E6ewquqTSd6Y1b1aK0nel+SsJOnu25LsT3JNkkNJfpTkxuHWNyR5W5KHq+rgcO69Q8XAP6iqy4fZfzPJ78ybhwALAADY9rr7+jn9neSd65z/66y/Pyvd/bZTnYcACwAAmNzUZdpPFwIsAABgcif5suBtT4AFAABMqiODBQAAMI5OsiQBljLtAAAAI5HBAgAAJmcPFgAAwFgEWAAAAGMoRS4AAABGsyQZLEUuAAAARiKDBQAATKu9BwsAAGA8S7JEUIAFAABsguXIYNmDBQAAMBIZLAAAYHqWCAIAAIxEgAUAADCCTqKKIAAAwDh6STJYilwAAACMRAYLAACY3pJksARYAADA9OzBAgAAGEfJYAEAAIygszRLBBW5AAAAGIkMFgAAMLGyBwsAAGA0S7JEUIAFAABMb0kCLHuwAAAARiKDBQAATG9JMlgCLAAAYFodRS4AAADGsiwvGh5lD1ZV7a2qZ6vqkQ36q6o+XFWHquqhqnr9GOMCAADbRC/Y5lgkJqmqq6rqiaHvlpnzL62qe6rqG8PnufPmMVaRi48nueoE/Vcn2TW0PUluHWlcAACA5AXGJFW1I8lHh/7LklxfVZcN99yS5N7u3pXk3uH4hEYJsLr7viTfO8El1ya5s1d9Kck5VXXhGGMDAAAsEJNckeRQdz/Z3T9O8qnh2uP33DF8vyPJm+fNY7PKtF+U5KmZ45XhHAAAsASqF2sj2CgmOVGsckF3H06S4fNl8wbZrCIX65UMWfcfU1XtyWrKLq94xSumnBMAALBZFq8ieF5VHZg5vr27bz+F+zeKSU46VjkZmxVgrSS5eOZ4Z5Kn17tw+Id0e5Ls3r17SWqNAADAGewkC1XM8Vx3717g/o1ikrM3OJ8kz1TVhd19eFhO+Oy8QTZrieC+JG8fKndcmeT7x1NtAAAAm2CjmOT+JLuq6tKqOjvJdcO1x++5Yfh+Q5LPzRtklAxWVX0yyRuzmrZbSfK+JGclSXfflmR/kmuSHEryoyQ3jjEuAACwTUy8Nu2FxiTdfaSqbk5yd5IdSfZ296PDY9+f5NNVdVOSbyV567x5jBJgdff1c/o7yTvHGAsAANh+pn7R8CIxSXfvz2oAtvb8d5O86VTmsVl7sAAAgGW2JNUVNmsPFgAAwBlPBgsAAJjekmSwBFgAAMCkRnxZ8GlPgAUAAExv8RcNbwsCLAAAYHpLksFS5AIAAGAkMlgAAMDk7MECAAAYiwALAABgBEtURdAeLAAAgJHIYAEAANNbkgyWAAsAAJieAAsAAGAc9mABAABwSgRYAAAAI7FEEAAAmN6SLBEUYAEAANNaovdgCbAAAIDpCbAAAABGsiQBliIXAAAAI5HBAgAAJlWxBwsAAGA8AiwAAIARLFEVQXuwAAAARiKDBQAATG9JMlgCLAAAYHoCLAAAgHEsyx4sARYAADC9JQmwFLkAAAAYiQwWAAAwrc7SZLAEWAAAwOSWZQ+WJYIAAMD0esE2R1VdVVVPVNWhqrplnf5zq+quqnqoqr5SVa8Zzv98VR2cac9X1buHvt+vqm/P9F0zbx4yWAAAwOSmzGBV1Y4kH03yq0lWktxfVfu6+7GZy96b5GB3v6WqXjVc/6bufiLJ5TPP+XaSu2bu+6Pu/sOTnYsMFgAAsN1dkeRQdz/Z3T9O8qkk16655rIk9yZJd389ySVVdcGaa96U5G+6++9e6EQEWAAAwPQWXyJ4XlUdmGl7Zp5+UZKnZo5XhnOzvpbkN5Okqq5I8sokO9dcc12ST645d/OwrHBvVZ07788UYAEAANNaNLhaDbCe6+7dM+32mRFqg1FnvT/JuVV1MMm/TPLVJEd+8oCqs5P8RpL/beaeW5P8XFaXEB5O8oF5f6o9WAAAwKQq60dAI1pJcvHM8c4kT89e0N3PJ7kxSaqqkvzt0I67OsmD3f3MzD0/+V5Vf5zk8/MmIoMFAABsd/cn2VVVlw6ZqOuS7Ju9oKrOGfqS5LeT3DcEXcddnzXLA6vqwpnDtyR5ZN5EZLAAAIDpTVhFsLuPVNXNSe5OsiPJ3u5+tKreMfTfluTVSe6sqqNJHkty0/H7q+pns1qB8HfWPPoPquryYfbfXKf/HxFgAQAAk5v6RcPdvT/J/jXnbpv5/sUkuza490dJ/uk65992qvMQYAEAANObOMA6XQiwAACA6S1JgKXIBQAAwEhksAAAgGn19HuwThcCLAAAYHpLEmCNskSwqq6qqieq6lBV3bJO/xur6vtVdXBovzfGuAAAwPZQvVjbLhbOYFXVjiQfzWrd+JUk91fVvu5+bM2lf9Xdv77oeAAAwDa0jYKkRYyRwboiyaHufrK7f5zkU0muHeG5AAAA28oYAdZFSZ6aOV4Zzq31S1X1tar6s6r6hY0eVlV7qupAVR34zne+M8L0AACArbYsSwTHCLBqnXNr/xE8mOSV3f3aJP8+yZ9u9LDuvr27d3f37vPPP3+E6QEAAFuqR2jbxBgB1kqSi2eOdyZ5evaC7n6+u384fN+f5KyqOm+EsQEAgO1AgHXS7k+yq6ouraqzk1yXZN/sBVX18qqq4fsVw7jfHWFsAACA08bCVQS7+0hV3Zzk7iQ7kuzt7ker6h1D/21JfivJ71bVkSR/n+S67t5GcSgAAPBCVbbXPqpFjPKi4WHZ3/41526b+f6RJB8ZYywAAGAbEmABAACMo5ZkAZsACwAAmNY2K1SxiDGKXAAAABAZLAAAYBMocgEAADAWARYAAMA4ZLAAAADGsiQBliIXAAAAI5HBAgAAptWWCAIAAIxHgAUAALC4yvJksOzBAgAAGIkMFgAAML1ejhSWAAsAAJjcsiwRFGABAADT6ihyAQAAMJY6ttUz2ByKXAAAAIxEBgsAAJjekiwRlMECAAAmV71Ym/v8qquq6omqOlRVt6zTf25V3VVVD1XVV6rqNTN936yqh6vqYFUdmDn/0qq6p6q+MXyeO28eAiwAAGBandUy7Yu0E6iqHUk+muTqJJclub6qLltz2XuTHOzuf5bk7Uk+tKb/l7v78u7ePXPuliT3dveuJPcOxyckwAIAACY3cQbriiSHuvvJ7v5xkk8luXbNNZdlNUhKd389ySVVdcGc516b5I7h+x1J3jxvIgIsAABgOzivqg7MtD0zfRcleWrmeGU4N+trSX4zSarqiiSvTLJz6Oskf1FVD6x57gXdfThJhs+XzZukIhcAAMD0Fi9y8dya5Xuz6iRGfH+SD1XVwSQPJ/lqkiND3xu6++mqelmSe6rq69193wuZpAALAACYVOXkClUsYCXJxTPHO5M8PXtBdz+f5MYkqapK8rdDS3c/PXw+W1V3ZXXJ4X1JnqmqC7v7cFVdmOTZeROxRBAAAJjWogUu5hS5SHJ/kl1VdWlVnZ3kuiT7Zi+oqnOGviT57ST3dffzVfWiqnrJcM2LkvxakkeG6/YluWH4fkOSz82biAwWAACwrXX3kaq6OcndSXYk2dvdj1bVO4b+25K8OsmdVXU0yWNJbhpuvyDJXatJrfx0kj/p7j8f+t6f5NNVdVOSbyV567y5CLAAAIDJTbxEMN29P8n+Nedum/n+xSS71rnvySSv3eCZ303yplOZhwALAACY3sQB1ulCgAUAAExu6gzW6UKABQAATKuTHFuOCEsVQQAAgJHIYAEAANNbjgSWAAsAAJiePVgAAABjmf+y4DOCAAsAAJjcsmSwFLkAAAAYiQwWAAAwrY4iFwAAAGOoJGUPFgAAwEiObfUENoc9WAAAACORwQIAACZniSAAAMAYFLkAAAAYS3vRMAAAwFi8aBgAAIBTMkqAVVVXVdUTVXWoqm5Zp7+q6sND/0NV9foxxgUAALaJ7sXaNrHwEsGq2pHko0l+NclKkvural93PzZz2dVJdg3tF5PcOnwCAABnuk7Ke7BO2hVJDnX3k9394ySfSnLtmmuuTXJnr/pSknOq6sIRxgYAALaDJclgjRFgXZTkqZnjleHcqV6TJKmqPVV1oKoOfOc73xlhegAAAJtjjACr1jm3NsQ8mWtWT3bf3t27u3v3+eefv/DkAACA00Av2LaJMcq0ryS5eOZ4Z5KnX8A1AADAGaq20TK/RYyRwbo/ya6qurSqzk5yXZJ9a67Zl+TtQzXBK5N8v7sPjzA2AACwHSzJHqyFM1jdfaSqbk5yd5IdSfZ296NV9Y6h/7Yk+5Nck+RQkh8luXHRcQEAgG2ikyxJFcExlgimu/dnNYiaPXfbzPdO8s4xxgIAADhdjRJgAQAAbKTSS7MHS4AFAABMT4AFAAAwEgEWAADACJaoyMUYZdoBAACIAAsAANgE1b1Qm/v8qquq6omqOlRVt6zTf25V3VVVD1XVV6rqNcP5i6vqC1X1eFU9WlXvmrnn96vq21V1cGjXzJuHJYIAAMD0JtyDVVU7knw0ya8mWUlyf1Xt6+7HZi57b5KD3f2WqnrVcP2bkhxJ8p7ufrCqXpLkgaq6Z+beP+ruPzzZuchgAQAAE+vVAGuRdmJXJDnU3U9294+TfCrJtWuuuSzJvUnS3V9PcklVXdDdh7v7weH8D5I8nuSiF/qXCrAAAIDt4LyqOjDT9sz0XZTkqZnjlfzjIOlrSX4zSarqiiSvTLJz9oKquiTJ65J8eeb0zcOywr1Vde68SQqwAACAaXXGyGA91927Z9rtMyPUBqPOen+Sc6vqYJJ/meSrWV0euPqAqhcn+UySd3f388PpW5P8XJLLkxxO8oF5f6o9WAAAwPSmLdO+kuTimeOdSZ6evWAImm5MkqqqJH87tFTVWVkNrj7R3Z+dueeZ49+r6o+TfH7eRGSwAACAyU1cRfD+JLuq6tKqOjvJdUn2/RfjV50z9CXJbye5r7ufH4KtjyV5vLs/uOaeC2cO35LkkXkTkcECAACmN2EVwe4+UlU3J7k7yY4ke7v70ap6x9B/W5JXJ7mzqo4meSzJTcPtb0jytiQPD8sHk+S93b0/yR9U1eVZXW74zSS/M28uAiwAAGDbGwKi/WvO3Tbz/YtJdq1z319n/T1c6e63neo8BFgAAMC0Osmx6TJYpxMBFgAAMLGTepfVGUGABQAATE+ABQAAMJIlCbCUaQcAABiJDBYAADAtRS4AAADG0kkf2+pJbAoBFgAAMD17sAAAADgVMlgAAMC07MECAAAY0ZIsERRgAQAA0xNgAQAAjKGXJsBS5AIAAGAkMlgAAMC0Oskx78ECAAAYx5IsERRgAQAA0xNgAQAAjKGX5j1YilwAAACMRAYLAACYVifdilwAAACMY0mWCAqwAACA6S1JkQt7sAAAAEYigwUAAEyr24uGAQAARrMkSwQFWAAAwORaBgsAAGAMvTQZLEUuAAAARiKDBQAATKvjPVgAAACjaXuwAAAAFtZJekkyWAvtwaqql1bVPVX1jeHz3A2u+2ZVPVxVB6vqwCJjAgAA20z3agZrkTZHVV1VVU9U1aGqumWd/nOr6q6qeqiqvlJVr5l378nGO7MWLXJxS5J7u3tXknuH4438cndf3t27FxwTAADgJ6pqR5KPJrk6yWVJrq+qy9Zc9t4kB7v7nyV5e5IPncS9pxLvJFk8wLo2yR3D9zuSvHnB5wEAAGegPtYLtTmuSHKou5/s7h8n+VRWY5VZl2U1SEp3fz3JJVV1wZx7TzneWTTAuqC7Dw+TPJzkZRtc10n+oqoeqKo9C44JAABsN9MuEbwoyVMzxyvDuVlfS/KbSVJVVyR5ZZKdc+492XjnJ+YWuaiqv0zy8nW6/sd59854Q3c/XVUvS3JPVX29u+/bYLw9SY4HYT+sqidOYRzmOy/Jc1s9CU6J32x78XttP36z7eeM/M3q3f96q6cwlTPy9zoNvHKrJ3AqfpD/++6/7P/9vAUf8zNr6jnc3t23D99rnevXpr3en+RDVXUwycNJvprkyEnee9LmBljd/Ssb9VXVM1V1YXcfrqoLkzy7wTOeHj6fraq7spqGWzfAGv4h3b5eH4urqgP2wW0vfrPtxe+1/fjNth+/2fbi9yJJuvuqiYdYSXLxzPHOJE+vmcPzSW5MkqqqJH87tJ89wb0nFe/MWnSJ4L4kNwzfb0jyubUXVNWLquolx78n+bUkjyw4LgAAwHH3J9lVVZdW1dlJrstqrPITVXXO0Jckv53kviHoOtG9c+OdtRZ9D9b7k3y6qm5K8q0kbx0m/18n+Q/dfU2SC5LctRok5qeT/El3//mC4wIAACRJur87IckAAAMNSURBVPtIVd2c5O4kO5Ls7e5Hq+odQ/9tSV6d5M6qOprksSQ3neje4dHrxjsnUt3L8cIvVlXVnpm1qmwDfrPtxe+1/fjNth+/2fbi92LZCLAAAABGsugeLAAAAAYCrCVUVW+tqker6lhVqepzmqqqq6rqiao6VFVz3xrO1qqqvVX1bFUp4rNNVNXFVfWFqnp8+O/Ed231nNhYVf1MVX2lqr42/F7/01bPiZNTVTuq6qtV9fmtngtsBgHWcnokqy9ZW7dUPluvqnYk+WiSq7P61vHrq+qyrZ0Vc3w8ydQlaBnXkSTv6e5XJ7kyyTv95+y09g9J/kV3vzbJ5Umuqqort3hOnJx3JXl8qycBm0WAtYS6+/Hu9gLn09sVSQ5195Pd/eMkn0py7RbPiRMYXp7+va2eByevuw9394PD9x9k9V8AL9raWbGRXvXD4fCsodlIfpqrqp1J/rsk/2Gr5wKbRYAFp6eLkjw1c7wS/+IHk6mqS5K8LsmXt3YmnMiw1OxgVl/0eU93+71Of/9Lkn+T5NhWTwQ2iwDrDFVVf1lVj6zTZEG2h1rnnP+nFiZQVS9O8pkk7x5eOMlpqruPdvflSXYmuaKqXrPVc2JjVfXrSZ7t7ge2ei6wmRZ90TCnqe7+la2eAwtZSXLxzPHOJE9v0VzgjFVVZ2U1uPpEd392q+fDyenu/6eq/q+s7ntUWOb09YYkv1FV1yT5mST/pKr+Y3f/91s8L5iUDBacnu5PsquqLq2qs5Ncl2TfFs8JzihVVUk+luTx7v7gVs+HE6uq86vqnOH7f5XkV5J8fWtnxYl097/t7p3dfUlW/3fs/xRcsQwEWEuoqt5SVStJfinJ/1FVd2/1nPgvdfeRJDcnuTurG+8/3d2Pbu2sOJGq+mSSLyb5+apaqaqbtnpOzPWGJG9L8i+q6uDQrtnqSbGhC5N8oaoeyur/CXVPdyv7DZx2qtu2DgAAgDHIYAEAAIxEgAUAADASARYAAMBIBFgAAAAjEWABAACMRIAFAAAwEgEWAADASARYAAAAI/n/AK0MR0vWVsgLAAAAAElFTkSuQmCC\n", "text/plain": [ "<Figure size 1152x432 with 2 Axes>" ] @@ -216,7 +216,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ @@ -237,7 +237,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 9, "metadata": {}, "outputs": [ { @@ -266,49 +266,63 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "domain_size: (2, 3)\n", - "direction: (-1, 1) , 5\n", - "break\n", - "write: [1, 0] read: [0, 2]\n", - "break\n", - "write: [1, 1] read: [0, 0]\n", - "write: [1, 2] read: [0, 1]\n", - "break\n", - "write: [0, 0] read: [1, 2]\n", - "break\n", - "write: [1, 0] read: [0, 2]\n", - "as tuples\n", - "((5, 33, 32),)\n", - "((5, 34, 30), (5, 35, 31))\n", - "((5, 30, 35),)\n", - "as lists\n", - "[[5, 33, 32]]\n", - "[[5, 34, 30], [5, 35, 31]]\n", - "[[5, 30, 35]]\n" + "domain_size: (4, 3)\n", + "[((1, 12, 14), (1, 15, 17), (1, 18, 20), (1, 21, 23)), ((2, 26, 24), (2, 29, 27), (2, 32, 30), (2, 35, 33)), ((3, 45, 36), (3, 46, 37), (3, 47, 38)), ((4, 48, 57), (4, 49, 58), (4, 50, 59)), ((5, 69, 62),), ((5, 70, 60), (5, 71, 61)), ((5, 60, 65), (5, 63, 68), (5, 66, 71)), ((5, 69, 62),), ((6, 72, 83),), ((6, 73, 81), (6, 74, 82)), ((6, 72, 83),), ((6, 75, 74), (6, 78, 77), (6, 81, 80)), ((7, 93, 85), (7, 94, 86)), ((7, 95, 84),), ((7, 86, 87), (7, 89, 90), (7, 92, 93)), ((7, 95, 84),), ((8, 96, 106), (8, 97, 107)), ((8, 98, 105),), ((8, 98, 105),), ((8, 101, 96), (8, 104, 99), (8, 107, 102))]\n", + "[[1, 12, 14], [1, 15, 17], [1, 18, 20], [1, 21, 23]]\n", + "[[2, 26, 24], [2, 29, 27], [2, 32, 30], [2, 35, 33]]\n", + "[[3, 45, 36], [3, 46, 37], [3, 47, 38]]\n", + "[[4, 48, 57], [4, 49, 58], [4, 50, 59]]\n", + "[[5, 69, 62]]\n", + "[[5, 70, 60], [5, 71, 61]]\n", + "[[5, 60, 65], [5, 63, 68], [5, 66, 71]]\n", + "[[6, 72, 83]]\n", + "[[6, 73, 81], [6, 74, 82]]\n", + "[[6, 75, 74], [6, 78, 77], [6, 81, 80]]\n", + "[[7, 93, 85], [7, 94, 86]]\n", + "[[7, 95, 84]]\n", + "[[7, 86, 87], [7, 89, 90], [7, 92, 93]]\n", + "[[8, 96, 106], [8, 97, 107]]\n", + "[[8, 98, 105]]\n", + "[[8, 101, 96], [8, 104, 99], [8, 107, 102]]\n" ] }, { "data": { "text/plain": [ - "[[[5, 33, 32]], [[5, 34, 30], [5, 35, 31]], [[5, 30, 35]]]" + "[[[1, 12, 14], [1, 15, 17], [1, 18, 20], [1, 21, 23]],\n", + " [[2, 26, 24], [2, 29, 27], [2, 32, 30], [2, 35, 33]],\n", + " [[3, 45, 36], [3, 46, 37], [3, 47, 38]],\n", + " [[4, 48, 57], [4, 49, 58], [4, 50, 59]],\n", + " [[5, 69, 62]],\n", + " [[5, 70, 60], [5, 71, 61]],\n", + " [[5, 60, 65], [5, 63, 68], [5, 66, 71]],\n", + " [[6, 72, 83]],\n", + " [[6, 73, 81], [6, 74, 82]],\n", + " [[6, 75, 74], [6, 78, 77], [6, 81, 80]],\n", + " [[7, 93, 85], [7, 94, 86]],\n", + " [[7, 95, 84]],\n", + " [[7, 86, 87], [7, 89, 90], [7, 92, 93]],\n", + " [[8, 96, 106], [8, 97, 107]],\n", + " [[8, 98, 105]],\n", + " [[8, 101, 96], [8, 104, 99], [8, 107, 102]]]" ] }, - "execution_count": 11, + "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "periodic_mapper = SparseLbPeriodicityMapper(method, mapping, domain_size, (True, True))\n", - "periodic_index_array = periodic_mapper.create_index_arr()\n", - "periodic_index_array" + "periodic_mapper.create_index_arr()\n", + "#periodic_index_array" ] }, { @@ -320,7 +334,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ @@ -362,9 +376,21 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "metadata": {}, - "outputs": [], + "outputs": [ + { + "ename": "NameError", + "evalue": "name 'index_arr' is not defined", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m<ipython-input-12-01d3d1672fc2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0;31m#index_field = ps.Field.create_from_numpy_array(\"idx\", index_arr, index_dimensions=1)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mindex_field\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mps\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mField\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcreate_generic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"idx\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mspatial_dimensions\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mindex_dimensions\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mindex_arr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mcollision_rule\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_collision_rule\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mNameError\u001b[0m: name 'index_arr' is not defined" + ] + } + ], "source": [ "#index_field = ps.Field.create_from_numpy_array(\"idx\", index_arr, index_dimensions=1)\n", "index_field = ps.Field.create_generic(\"idx\", spatial_dimensions=1, index_dimensions=1, dtype=index_arr.dtype)\n",