post-merge fixes; support node interpolation to fixed_mesh_model.py a…

…nd add it to the 3D simulation notebook;

docs/notebooks/3D_simulation.ipynb

@@ -220,6 +220,55 @@
     "model.simulator.run(save=True,purge=True)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "# interpolate and extrapolate the missing nodes\n",
+    "# find nearby nodes from the array indexer_full_sim, which is sorted by x-index\n",
+    "import itertools\n",
+    "from subsheat3D.fixed_mesh_model import interpolateNode\n",
+    "for ni in range(len(model.builder.nodes)):\n",
+    "    for nj in range(len(model.builder.nodes[ni])):\n",
+    "        if model.builder.nodes[ni][nj] is False:\n",
+    "            closest_x_up = []\n",
+    "            for j in range(ni,len(model.builder.nodes[nj])):\n",
+    "                matching_x = [ i[0] for i in model.builder.indexer_full_sim if i[0]==j ]\n",
+    "                closest_x_up = closest_x_up + list(set(matching_x))\n",
+    "                if len(matching_x)>0:\n",
+    "                    break\n",
+    "            closest_x_down = []\n",
+    "            for j in range(ni-1,-1,-1):\n",
+    "                matching_x = [ i[0] for i in model.builder.indexer_full_sim if i[0]==j ]\n",
+    "                closest_x_down = closest_x_down + list(set(matching_x))\n",
+    "                if len(matching_x)>0:\n",
+    "                    break\n",
+    "            closest_y_up = []\n",
+    "            for j in range(nj,len(model.builder.nodes[ni])):\n",
+    "                matching_y = [ i[1] for i in model.builder.indexer_full_sim if (i[1]==j and ((i[0] in closest_x_up) or i[0] in closest_x_down)) ]\n",
+    "                closest_y_up = closest_y_up + list(set(matching_y))\n",
+    "                if len(matching_y)>0:\n",
+    "                    break\n",
+    "            closest_y_down = []\n",
+    "            for j in range(nj-1,-1,-1):\n",
+    "                matching_y = [ i[1] for i in model.builder.indexer_full_sim if (i[1]==j and (i[0] in closest_x_up or i[0] in closest_x_down) ) ]\n",
+    "                closest_y_down = closest_y_down + list(set(matching_y))\n",
+    "                if len(matching_y)>0:\n",
+    "                    break\n",
+    "\n",
+    "            interpolationNodes = [  model.builder.nodes[i[0]][i[1]] for i in itertools.product(closest_x_up+closest_x_down, closest_y_up+closest_y_down)  ]\n",
+    "            interpolationNodes = [nn for nn in interpolationNodes if nn is not False]\n",
+    "            node = interpolateNode(interpolationNodes)\n",
+    "            node.X, node.Y = model.builder.grid.location_grid[ni,nj,:]\n",
+    "            model.builder.nodes[ni][nj] = node\n",
+    "        else:\n",
+    "            node = interpolateNode([model.builder.nodes[ni][nj]])  # \"interpolate\" the node from itself to make sure the same member variables exist at the end\n",
+    "            model.builder.nodes[ni][nj] = node\n"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 15,
@@ -246,7 +295,17 @@
    ],
    "source": [
     "from warmth.mesh_model import run\n",
-    "run(model,start_time=model.parameters.time_start,end_time=0)"
+    "import os\n",
+    "try:\n",
+    "    os.mkdir('mesh')\n",
+    "except FileExistsError:\n",
+    "    pass\n",
+    "try:\n",
+    "    os.mkdir('temp')\n",
+    "except FileExistsError:\n",
+    "    pass\n",
+    "run(model,start_time=model.parameters.time_start,end_time=0)\n",
+    "\n"
    ]
   },
   {

subsheat3D/fixed_mesh_model.py

@@ -1,7 +1,9 @@
 import numpy as np
 from mpi4py import MPI
 from scipy.interpolate import LinearNDInterpolator
+from typing import Iterator, List, Literal
 
+from warmth.build import single_node
 from warmth.logging import logger
 from subsheat3D.Helpers import  NodeParameters1D, top_crust,top_sed,thick_crust,  top_lith, top_asth, top_sed_id, bottom_sed_id
 from subsheat3D.resqpy_helpers import write_tetra_grid_with_properties, write_hexa_grid_with_properties
@@ -15,6 +17,46 @@
 #      NOTE: the default values of surface and base heat flow imply zero RHP in the crust
 #
 
+
+def interpolateNode(interpolationNodes: List[single_node], interpolationWeights=None) -> single_node:
+    assert len(interpolationNodes)>0
+    if interpolationWeights is None:
+        interpolationWeights = np.ones([len(interpolationNodes),1])
+    assert len(interpolationNodes)==len(interpolationWeights)
+    wsum = np.sum(np.array(interpolationWeights))
+    iWeightNorm = [ w/wsum for w in interpolationWeights]
+
+    node = single_node()
+    node.__dict__.update(interpolationNodes[0].__dict__)
+    node.X = np.sum( np.array( [node.X * w for node,w in zip(interpolationNodes,iWeightNorm)] ) ) 
+    node.Y = np.sum( np.array( [node.Y * w for node,w in zip(interpolationNodes,iWeightNorm)] ) )
+
+    times = range(node.result._depth.shape[1])
+    if node.subsidence is None:
+        node.subsidence = np.sum( np.array( [ [node.result.seabed(t) for t in times] * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) 
+    if node.crust_ls is None:
+        node.crust_ls = np.sum( np.array( [ [node.result.crust_thickness(t) for t in times] * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) 
+    if node.lith_ls is None:
+        node.lith_ls = np.sum( np.array( [ [node.result.lithosphere_thickness(t) for t in times] * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) 
+
+    if node.beta is None:
+        node.beta = np.sum( np.array( [node.beta * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) 
+    if node.kAsth is None:
+        node.kAsth = np.sum( np.array( [node.kAsth * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) 
+    if node.kLith is None:
+        node.kLith = np.sum( np.array( [node.kLith * w for node,w in zip(interpolationNodes,iWeightNorm)] ) , axis = 0) 
+    if node._depth_out is None:
+        node._depth_out = np.sum([node.result._depth_out*w for n,w in zip(interpolationNodes[0:1], [1] )], axis=0)
+    if node.temperature_out is None:
+        node.temperature_out = np.sum([n.result.temperature_out*w for n,w in zip(interpolationNodes[0:1], [1] )], axis=0)
+
+    if node.sed is None:
+        node.sed = np.sum([n.sed*w for n,w in zip(interpolationNodes,iWeightNorm)], axis=0)
+    if node.sed_thickness_ls is None:
+        node.sed_thickness_ls =  node.sed[-1,1,:] - node.sed[0,0,:]    
+    return node
+
+
 #
 class UniformNodeGridFixedSizeMeshModel:
     """Manages a 3D heat equation computation using dolfinx

warmth/mesh_model.py

@@ -10,7 +10,7 @@
 from .model import Model
 from warmth.logging import logger
 from .mesh_utils import  NodeParameters1D, top_crust,top_sed,thick_crust,  top_lith, top_asth, top_sed_id, bottom_sed_id,NodeGrid
-from .resqpy_helpers import write_tetra_grid_with_properties, write_hexa_grid_with_properties,read_mesh_resqml
+from .resqpy_helpers import write_tetra_grid_with_properties, write_hexa_grid_with_properties,read_mesh_resqml_hexa
 def tic():
     #Homemade version of matlab tic and toc functions
     import time
@@ -1374,4 +1374,4 @@ def run( model:Model, run_simulation=True, start_time=182, end_time=0, out_dir =
     print("total time solve: " , time_solve)
     EPCfilename = mm2.write_hexa_mesh_resqml("temp/")
     print("RESQML model written to: " , EPCfilename)
-    read_mesh_resqml(EPCfilename)  # test reading of the .epc file
+    read_mesh_resqml_hexa(EPCfilename)  # test reading of the .epc file