diff --git a/examples/example5/example5_withSOCE.ipynb b/examples/example5/example5_withSOCE.ipynb
new file mode 100644
index 0000000..c70075b
--- /dev/null
+++ b/examples/example5/example5_withSOCE.ipynb
@@ -0,0 +1,838 @@
+{
+ "cells": [
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "f65f18d7",
+ "metadata": {},
+ "source": [
+ "# Example 5: Generic cell signaling system in 3D, illustrating a SOCE-like effect\n",
+ "\n",
+ "This example is slightly altered from the normal version of example 5 to illustrate the role of store-operated calcium entry (SOCE) accompanying calcium release from an interior compartment.\n",
+ "\n",
+ "Geometry is divided into 4 domains; two volumes, and two surfaces:\n",
+ "- cytosol (Cyto): $\\Omega_{Cyto}$\n",
+ "- endoplasmic reticulum volume (ER): $\\Omega_{ER}$\n",
+ "- plasma membrane (PM): $\\Gamma_{PM}$\n",
+ "- ER membrane (ERm): $\\Gamma_{ERm}$\n",
+ "\n",
+ "For simplicity, here we consider a \"cube-within-a-cube\" geometry, in which the smaller\n",
+ "inner cube represents a section of ER and one face of the outer cube ($x=0$) represents the PM. The other\n",
+ "faces of the outer cube are treated as no flux boundaries. The space outside\n",
+ "the inner cube but inside the outer cube is classified as cytosol.\n",
+ "\n",
+ "There are three function-spaces on these three domains:\n",
+ "\n",
+ "$$\n",
+ "u^{Cyto} = [A, B] \\quad \\text{on} \\quad \\Omega^{Cyto}\\\\\n",
+ "u^{ER} = [AER] \\quad \\text{on} \\quad \\Omega^{ER}\\\\\n",
+ "v^{ERm} = [R, Ro] \\quad \\text{on} \\quad \\Gamma^{ERm}\n",
+ "$$\n",
+ "\n",
+ "In words, this says that species A and B reside in the cytosolic volume, \n",
+ "species AER corresponds to an amount of species A that lives in the ER volume,\n",
+ "and species R (closed receptor/channel) and Ro (open receptor/channel) reside on the ER membrane.\n",
+ "\n",
+ "In this model, species B reacts with a receptor/channel, R, on the ER membrane, causing it to open (change state from R->Ro), \n",
+ "allowing species A to flow out of the ER and into the cytosol. \n",
+ "Note that this is roughly similar to an IP3 pulse at the PM, leading to Ca2+ release from the ER,\n",
+ "where, by analogy, species B is similar to IP3 and species A is similar to Ca2+. A more comprehensive\n",
+ "model of Ca2+ dynamics in particular is implemented in Example 6.\n",
+ "\n",
+ "Here, we also introduce a flux analogous to SOCE, wherein the flux of A through the PM depends on depletion of A from the ER.\n",
+ "For simplicity, we scale this flux with the distance between the PM and ER ($d_{PM-ER}$) and the concentration of A at closest point of ER to the PM ($A_{ER,close}$):\n",
+ "\n",
+ "$$\n",
+ "J_{SOCE} = \\frac{d_0}{d_{PM-ER} + d_0} \\frac{1}{1 + \\left(\\frac{A_{ER}}{A_{ref}}\\right)^4},\n",
+ "$$\n",
+ "\n",
+ "where $d_0$ is a characteristic length scale controlling the likelihood of STIM-Orai binding and $A_{ref}$ governs the ER calcium depletion required to trigger SOCE."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a5edff5d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from matplotlib import pyplot as plt\n",
+ "import matplotlib.image as mpimg\n",
+ "img_A = mpimg.imread('example5-diagram.png')\n",
+ "plt.imshow(img_A)\n",
+ "plt.axis('off')"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "f543d75e",
+ "metadata": {},
+ "source": [
+ "As specified in our [mathematical documentation](https://rangamanilabucsd.github.io/smart/docs/math.html), assuming diffusive transport, the PDE and boundary condition for each of these volumetric species takes the form:\n",
+ "\n",
+ "$$\n",
+ "\\partial_t u_i^m - \\nabla \\cdot ( D_i^m \\nabla (u_i^m) ) - f_i^m(u_i^m) = 0 \\qquad \\text{ in } \\Omega^m\\\\\n",
+ "D_i \\nabla u_i^m \\cdot n^m - R_i^{q} (u^m, u^n, v^q) = 0 \\qquad \\text{ on } \\Gamma^{q}\n",
+ "$$\n",
+ "\n",
+ "and the surface species take the form:\n",
+ "\n",
+ "$$\n",
+ "\\partial_t v_i^q - \\nabla_S \\cdot (D_i^q \\nabla_S v ) - g_i^q ( u^m, u^n, v^q ) = 0 \\qquad \\text{ on } \\Gamma^{q}\\\\\n",
+ "D_i \\nabla v_i^q \\cdot n^q = 0 \\qquad \\text{ on } \\partial\\Gamma^{q}\n",
+ "$$\n",
+ "\n",
+ "Our reaction terms and boundary conditions are chosen according to the system described above. For the purposes of this simplified example we use linear mass action in all reaction terms except SOCE. Explicitly writing out this system of PDEs, we have:\n",
+ "\n",
+ "\\begin{align}\n",
+ " \\partial_t u_B^{Cyto} - D_B^{Cyto} \\nabla^2 u_B^{Cyto} + k_{2f} u_B^{Cyto} &= 0 \\qquad \\text{ in } \\Omega^{Cyto}\\\\\n",
+ " D_B^{Cyto} \\nabla u_B^{Cyto} \\cdot n^{Cyto} + j_1[t] &= 0 \\qquad \\text{ on } \\Gamma^{PM} \\nonumber\\\\\n",
+ " D_B^{Cyto} \\nabla u_B^{Cyto} \\cdot n^{Cyto} + k_{3f} v_R^{ERm} u_B^{Cyto} - k_{3r} v_{Ro}^{ERm} &= 0 \\qquad \\text{ on } \\Gamma^{ERm} \\nonumber\\\\\n",
+ " \\nonumber \\\\\n",
+ " \\partial_t u_A^{Cyto} - D_A^{Cyto} \\nabla^2 u_A^{Cyto} &= 0 \\qquad \\text{ in } \\Omega^{Cyto}\\\\\n",
+ " D_A^{Cyto} \\nabla u_A^{Cyto} \\cdot n^{Cyto} - J_{SOCE} &= 0 \\qquad \\text{ on } \\Gamma^{PM} \\nonumber\\\\\n",
+ " D_A^{Cyto} \\nabla u_A^{Cyto} \\cdot n^{Cyto} - k_{4,Vmax} v_{Ro}^{ERm} (u_{AER}^{ER} - u_A^{Cyto}) &= 0 \\qquad \\text{ on } \\Gamma^{ERm} \\nonumber\\\\\n",
+ " \\nonumber \\\\\n",
+ " \\partial_t u_{AER}^{ER} - D_{AER}^{ER} \\nabla^2 u_{AER}^{ER} &= 0 \\qquad \\text{ in } \\Omega^{ER}\\\\\n",
+ " D_{AER}^{ER} \\nabla u_{AER}^{ER} \\cdot n^{ER} + k_{4,Vmax} v_{Ro}^{ERm} (u_{AER}^{ER} - u_A^{Cyto}) &= 0 \\qquad \\text{ on } \\Gamma^{ERm} \\nonumber\\\\\n",
+ " \\nonumber \\\\\n",
+ " \\partial_t v_{R}^{ERm} - D_{R}^{ERm} \\nabla^2 v_{R}^{ERm} - \n",
+ " k_{3f} v_R^{ERm} u_B^{Cyto} + k_{3r} v_{Ro}^{ERm} &= 0 \\qquad \\text{ on } \\Gamma^{ERm}\\\\\n",
+ " \\nonumber \\\\\n",
+ " \\partial_t v_{Ro}^{ERm} - D_{Ro}^{ERm} \\nabla^2 v_{Ro}^{ERm} +\n",
+ " k_{3f} v_R^{ERm} u_B^{Cyto} - k_{3r} v_{Ro}^{ERm} &= 0 \\qquad \\text{ on } \\Gamma^{ERm}\\\\\n",
+ "\\end{align}"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "28f77cbf",
+ "metadata": {},
+ "source": [
+ "## Code imports and initialization"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "cc398816",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import logging\n",
+ "\n",
+ "import dolfin as d\n",
+ "import sympy as sym\n",
+ "import numpy as np\n",
+ "import pathlib\n",
+ "import gmsh # must be imported before pyvista if dolfin is imported first\n",
+ "\n",
+ "from smart import config, mesh, model, mesh_tools, visualization\n",
+ "from smart.model_assembly import (\n",
+ " Compartment,\n",
+ " Parameter,\n",
+ " Reaction,\n",
+ " Species,\n",
+ " SpeciesContainer,\n",
+ " ParameterContainer,\n",
+ " CompartmentContainer,\n",
+ " ReactionContainer,\n",
+ ")\n",
+ "from smart.units import unit"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "c8650536",
+ "metadata": {},
+ "source": [
+ "We will set the logging level to `INFO`. This will display some output during the simulation. If you want to get even more output you could set the logging level to `DEBUG`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9ed0899a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "logger = logging.getLogger(\"smart\")\n",
+ "logger.setLevel(logging.INFO)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "defc1095",
+ "metadata": {},
+ "source": [
+ "Futhermore, you could also save the logs to a file by attaching a file handler to the logger as follows.\n",
+ "\n",
+ "```\n",
+ "file_handler = logging.FileHandler(\"filename.log\")\n",
+ "file_handler.setFormatter(logging.Formatter(smart.config.base_format))\n",
+ "logger.addHandler(file_handler)\n",
+ "```"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "95b9d865",
+ "metadata": {},
+ "source": [
+ "First, we define the various units for the inputs"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4f4023cf",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Aliases - base units\n",
+ "uM = unit.uM\n",
+ "um = unit.um\n",
+ "molecule = unit.molecule\n",
+ "sec = unit.sec\n",
+ "# Aliases - units used in model\n",
+ "D_unit = um**2 / sec\n",
+ "flux_unit = molecule / (um**2 * sec)\n",
+ "vol_unit = uM\n",
+ "surf_unit = molecule / um**2"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "46582d26",
+ "metadata": {},
+ "source": [
+ "## Generate model\n",
+ "Next we generate the model described in the equations above.\n",
+ "\n",
+ "### Compartments\n",
+ "As described above, the three compartments are the cytosol (\"Cyto\"), the plasma membrane (\"PM\"), the ER membrane (\"ERm\"), and the ER interior volume (\"ER\").\n",
+ "\n",
+ "Note that, as shown, we can also specify nonadjacency for compartments; this is not strictly necessary, but will generally speed up the simulations."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b2f34b3c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Cyto = Compartment(\"Cyto\", 3, um, 1)\n",
+ "PM = Compartment(\"PM\", 2, um, 10)\n",
+ "ER = Compartment(\"ER\", 3, um, 2)\n",
+ "ERm = Compartment(\"ERm\", 2, um, 12)\n",
+ "PM.specify_nonadjacency(['ERm', 'ER'])\n",
+ "ERm.specify_nonadjacency(['PM'])"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "d93ce862",
+ "metadata": {},
+ "source": [
+ "Initialize a compartment container and add the 4 compartments to it."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "701577e8",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "cc = CompartmentContainer()\n",
+ "cc.add([ERm, ER, PM, Cyto])"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "0a6acf0b",
+ "metadata": {},
+ "source": [
+ "### Species\n",
+ "In this case, we have 5 species across 3 different compartments. We define each in turn:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8d991278",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "A = Species(\"A\", 0.01, vol_unit, 1.0, D_unit, \"Cyto\")\n",
+ "B = Species(\"B\", 0.0, vol_unit, 1.0, D_unit, \"Cyto\")\n",
+ "AER = Species(\"AER\", 200.0, vol_unit, 5.0, D_unit, \"ER\")\n",
+ "# Uniform initial condition of R for simplicity\n",
+ "R1 = Species(\"R1\", 1.0, surf_unit, 0.02, D_unit, \"ERm\")\n",
+ "R1o = Species(\"R1o\", 0.0, surf_unit, 0.02, D_unit, \"ERm\")"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "b60826cc",
+ "metadata": {},
+ "source": [
+ "Create species container and add the 5 species objects to it."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e86bebaf",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "sc = SpeciesContainer()\n",
+ "sc.add([R1o, R1, AER, B, A])"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "af900a73",
+ "metadata": {},
+ "source": [
+ "### Parameters and Reactions\n",
+ "\n",
+ "Parameters and reactions are generally defined together, although the order does not strictly matter. We define them in turn as follows:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d2120ac6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Degradation of B in the cytosol\n",
+ "k2f = Parameter(\"k2f\", 10, 1 / sec)\n",
+ "r2 = Reaction(\n",
+ " \"r2\", [\"B\"], [], param_map={\"on\": \"k2f\"}, reaction_type=\"mass_action_forward\"\n",
+ ")\n",
+ "\n",
+ "# Activating receptors on ERm with B\n",
+ "k3f = Parameter(\"k3f\", 100, 1 / (uM * sec))\n",
+ "k3r = Parameter(\"k3r\", 100, 1 / sec)\n",
+ "r3 = Reaction(\"r3\", [\"B\", \"R1\"], [\"R1o\"], {\"on\": \"k3f\", \"off\": \"k3r\"})\n",
+ "# Release of A from ERm to cytosol\n",
+ "k4Vmax = Parameter(\"k4Vmax\", 2000, 1 / (uM * sec))\n",
+ "r4 = Reaction(\n",
+ " \"r4\",\n",
+ " [\"AER\"],\n",
+ " [\"A\"],\n",
+ " param_map={\"Vmax\": \"k4Vmax\"},\n",
+ " species_map={\"R1o\": \"R1o\", \"uER\": \"AER\", \"u\": \"A\"},\n",
+ " eqn_f_str=\"Vmax*R1o*(uER-u)\",\n",
+ ")"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "7fdd58b8",
+ "metadata": {},
+ "source": [
+ "We define an additional reaction as the time-dependent production of species B at the plasma membrane. In this case, we define a pulse-type function as the derivative of an arctan function. Note that this is useful because we can provide an expression to use for pre-integration.\n",
+ "\n",
+ "$$\n",
+ "j_{int}[t] = V_{max} \\arctan\\left({m (t - t_0)}\\right)\\\\\n",
+ "j_1[t] = \\frac{dj_{int}[t]}{dt} = \\frac{m V_{max}}{1 + m^2 (t-t_0)^2}\n",
+ "$$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "51abf270",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Vmax, t0, m = 500, 0.1, 200\n",
+ "t = sym.symbols(\"t\")\n",
+ "pulseI = Vmax * sym.atan(m * (t - t0))\n",
+ "pulse = sym.diff(pulseI, t)\n",
+ "j1pulse = Parameter.from_expression(\n",
+ " \"j1pulse\", pulse, flux_unit, use_preintegration=True, preint_sym_expr=pulseI\n",
+ ")\n",
+ "r1 = Reaction(\n",
+ " \"r1\",\n",
+ " [],\n",
+ " [\"B\"],\n",
+ " param_map={\"J\": \"j1pulse\"},\n",
+ " eqn_f_str=\"J\",\n",
+ " explicit_restriction_to_domain=\"PM\",\n",
+ ")"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "b5e218b0",
+ "metadata": {},
+ "source": [
+ "We can plot the time-dependent input by converting the sympy expression to a numpy function using lambdify."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a0829a6b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from sympy.utilities.lambdify import lambdify\n",
+ "pulse_func = lambdify(t, pulse, 'numpy') # returns a numpy-ready function\n",
+ "tArray = np.linspace(0, 1, 100)\n",
+ "pulse_vals = pulse_func(tArray)\n",
+ "plt.plot(tArray, pulse_vals)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "55c27105",
+ "metadata": {},
+ "source": [
+ "Finally, we define the store-operated calcium (or here, A) influx. \n",
+ "This is not immediately straightforward, as a flux at the PM depends on concentration of A in the ER.\n",
+ "One option would be to explicitly introduce contacts between the ER and PM, but this would imply direct entry of calcium into the ER through Orai1.\n",
+ "Instead, we assume the likelihood of STIM-Orai1 binding scales with PM-ER distance (if they are very close, membrane fluctuations are likely to allow contact between the ERM and PM occasionally, for instance).\n",
+ "To adopt this strategy, we must construct a map between each point on the PM and the closest point on the ER, which we do after initializing the model below.\n",
+ "For now, we define the flux as a parameter varying over the PM mesh by calling `Parameter.mesh_quantity()`\n",
+ ""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a56bbccb",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "PSOCE = Parameter.mesh_quantity(\"PSOCE\", 0, unit.dimensionless, compartment=\"PM\")\n",
+ "J0_SOCE = Parameter(\"J0_SOCE\", 1e4, flux_unit)\n",
+ "rSOCE = Reaction(\n",
+ " \"rSOCE\",\n",
+ " [],\n",
+ " [\"A\"],\n",
+ " param_map={\"P\": \"PSOCE\", \"J0\": \"J0_SOCE\"},\n",
+ " eqn_f_str=\"P*J0\",\n",
+ " explicit_restriction_to_domain=\"PM\",\n",
+ ")"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "3d883c6e",
+ "metadata": {},
+ "source": [
+ "Create containers for parameters and reactions and add all the parameters and reaction objects to them."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8ae2c2c1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pc = ParameterContainer()\n",
+ "rc = ReactionContainer()\n",
+ "pc.add([k4Vmax, k3r, k3f, k2f, j1pulse, PSOCE, J0_SOCE])\n",
+ "rc.add([r1, r2, r3, r4, rSOCE])"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "15c35d39",
+ "metadata": {},
+ "source": [
+ "## Create/load in mesh\n",
+ "\n",
+ "In SMART we have different levels of meshes:\n",
+ "- Parent mesh: contains the entire geometry of the problem, including all surfaces and volumes\n",
+ "- Child meshes: submeshes (sections of the parent mesh) associated with individual compartments. Here, the child meshes are:\n",
+ " - Cyto: the portion of the outer cube outside of the inner cube, defined by `cell_markers = 1`\n",
+ " - ER: the inside portion of the inner cube, defined by `cell_markers = 2`\n",
+ " - PM: surface mesh where x=0, defined by `facet_markers = 10`\n",
+ " - ERm: surface mesh corresponding to all faces of the inner cube, defined by `facet_markers = 12`\n",
+ "\n",
+ "Here we create a UnitCube mesh as the Parent mesh, defined by\n",
+ "\n",
+ "$$\n",
+ "\\Omega = [0, 1] \\times [0, 1] \\times [0, 1] \\subset \\mathbb{R}^3\n",
+ "$$\n",
+ "\n",
+ "\n",
+ "The ER is a cube within the exterior cube, with dimensions 0.4 by 0.4 by 0.4 and a tunable gap between the ER and PM to test for different strengths of SOCE below."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e77f97fe",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ER_PM_gap = 0.05 # must be greater than 0 and less than 0.6\n",
+ "ER_PM_gap = 0.05 * round(ER_PM_gap/0.05) # round to ensure flat edges in mesh\n",
+ "def cur_cube_condition(cell, xmin=ER_PM_gap, xmax=ER_PM_gap+0.4):\n",
+ " \"\"\"\n",
+ " Returns true when inside an inner cube region defined as:\n",
+ " xmin <= x <= xmax, ymin <= y <= ymax, zmin <= z <= zmax\n",
+ " \"\"\"\n",
+ " ymin = 0.3\n",
+ " ymax = 0.7\n",
+ " zmin = 0.3\n",
+ " zmax = 0.7\n",
+ " return (\n",
+ " (xmin - d.DOLFIN_EPS < cell.midpoint().x() < xmax + d.DOLFIN_EPS)\n",
+ " and (ymin - d.DOLFIN_EPS < cell.midpoint().y() < ymax + d.DOLFIN_EPS)\n",
+ " and (zmin - d.DOLFIN_EPS < cell.midpoint().z() < zmax + d.DOLFIN_EPS)\n",
+ " )\n",
+ "domain, facet_markers, cell_markers = mesh_tools.create_cubes(condition=cur_cube_condition, N=20)\n",
+ "visualization.plot_dolfin_mesh(domain, cell_markers, clip_plane=(1,\n",
+ " 1, 0), clip_origin=(0.5, 0.5, 0.5))\n"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "035113ed",
+ "metadata": {},
+ "source": [
+ "By default, `smart.mesh_tools.create_cubes` marks all faces of the outer cube as \"10\", our marker value associated with PM. Here, since we are only treating the x=0 face as PM, we alter the facet markers on all other faces, setting them equal to zero. They are then treated as no-flux boundaries not belonging to a designated surface compartment."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "fe56e162",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for face in d.faces(domain):\n",
+ " if face.midpoint().x() > d.DOLFIN_EPS and facet_markers[face] == 10:\n",
+ " facet_markers[face] = 0"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "c17ebcbb",
+ "metadata": {},
+ "source": [
+ "We now save the mesh as an h5 file and then read it into SMART as a `ParentMesh` object. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5657aab1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "mesh_folder = pathlib.Path(\"mesh\")\n",
+ "mesh_folder.mkdir(exist_ok=True)\n",
+ "mesh_path = mesh_folder / \"DemoCuboidsMesh.h5\"\n",
+ "mesh_tools.write_mesh(\n",
+ " domain, facet_markers, cell_markers, filename=mesh_path\n",
+ ")\n",
+ "parent_mesh = mesh.ParentMesh(\n",
+ " mesh_filename=str(mesh_path),\n",
+ " mesh_filetype=\"hdf5\",\n",
+ " name=\"parent_mesh\",\n",
+ ")"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "0f1cf8a4",
+ "metadata": {},
+ "source": [
+ "## Model and solver initialization\n",
+ "\n",
+ "Now we are ready to set up the model. First we load the default configurations and set some configurations for the current solver."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "8c1a2a92",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "conf = config.Config()\n",
+ "conf.solver.update(\n",
+ " {\n",
+ " \"final_t\": 1,\n",
+ " \"initial_dt\": 0.01,\n",
+ " \"time_precision\": 6,\n",
+ " }\n",
+ ")"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "1511acc7",
+ "metadata": {},
+ "source": [
+ "We create a model using the different containers and the parent mesh. For later reference, we save the model information as a pickle file. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e3c3c30f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "model_cur = model.Model(pc, sc, cc, rc, conf, parent_mesh)\n",
+ "model_cur.to_pickle('model_cur.pkl')"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "dfb70c74",
+ "metadata": {},
+ "source": [
+ "Note that we could later load the model information from the pickle file using the line:\n",
+ "```Python\n",
+ "model_cur = model.from_pickle(model_cur.pkl)\n",
+ "```\n",
+ "\n",
+ "Next we need to initialize the model and solver."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c8976aa2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "model_cur.initialize()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "70f79003",
+ "metadata": {},
+ "source": [
+ "Now create mapping between PM and closest point on the ER surface for spatially dependent SOCE."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "db8d9dca",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "Pfun = PSOCE.dolfin_function # function associated with PSOCE values\n",
+ "V_PM = Pfun.function_space() # PM function space\n",
+ "# pull coordinate lists for PM and ERM\n",
+ "PMcoord = V_PM.tabulate_dof_coordinates()\n",
+ "ERMcoord = ERm.dolfin_mesh.coordinates()\n",
+ "# initialize vectors for distances and AER indices\n",
+ "AER_idx = np.zeros_like(Pfun.vector())\n",
+ "dist_vals = np.zeros_like(Pfun.vector())\n",
+ "# get parent mesh and vertex mapping to meshviews\n",
+ "mesh_ref = model_cur.parent_mesh.dolfin_mesh\n",
+ "ERM_map = ERm.dolfin_mesh.topology().mapping()[mesh_ref.id()].vertex_map()\n",
+ "ER_map = ER.dolfin_mesh.topology().mapping()[mesh_ref.id()].vertex_map()\n",
+ "# find closest ERM point for each PM coordinate\n",
+ "for i in range(len(PMcoord)):\n",
+ " curCoord = PMcoord[i]\n",
+ " dists = np.sqrt((curCoord[0]-ERMcoord[:,0])**2 + \n",
+ " (curCoord[1]-ERMcoord[:,1])**2 + \n",
+ " (curCoord[2]-ERMcoord[:,2])**2)\n",
+ " ERMidx = np.argmin(dists)\n",
+ " dist_vals[i] = dists[ERMidx]\n",
+ " global_idx = ERM_map[ERMidx]\n",
+ " ER_idx = np.nonzero(np.array(ER_map)==global_idx)[0]\n",
+ " if len(ER_idx) != 1:\n",
+ " raise ValueError(\"Node not found in ER mesh\")\n",
+ " AER_idx[i] = d.vertex_to_dof_map(AER.V)[ER_idx][0]\n",
+ "# now define PSOCE function values according to dist vals and AER values\n",
+ "d0 = 0.05\n",
+ "cref = 100.0\n",
+ "vals_new = (d0/(dist_vals+d0)) * (1/(1+(AER.u[\"u\"].vector()[AER_idx]/cref)**4))\n",
+ "PSOCE.dolfin_function.vector()[:] = vals_new\n",
+ "PSOCE.dolfin_function.vector().apply(\"insert\")"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "d05a1a75",
+ "metadata": {},
+ "source": [
+ "## Solving system and storing data\n",
+ "\n",
+ "We create some XDMF files where we will store the output "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5df9509e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Write initial condition(s) to file\n",
+ "results = dict()\n",
+ "result_folder = pathlib.Path(f\"results_ER_PM_gap_{ER_PM_gap}\")\n",
+ "result_folder.mkdir(exist_ok=True)\n",
+ "for species_name, species in model_cur.sc.items:\n",
+ " results[species_name] = d.XDMFFile(\n",
+ " model_cur.mpi_comm_world, str(result_folder / f\"{species_name}.xdmf\")\n",
+ " )\n",
+ " results[species_name].parameters[\"flush_output\"] = True\n",
+ " results[species_name].write(model_cur.sc[species_name].u[\"u\"], model_cur.t)\n",
+ "SOCEfile = d.XDMFFile(model_cur.mpi_comm_world, str(result_folder / f\"PSOCE.xdmf\"))\n",
+ "SOCEfile.parameters[\"flush_output\"] = True\n",
+ "SOCEfile.write(PSOCE.dolfin_function, model_cur.t)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "183376e6",
+ "metadata": {},
+ "source": [
+ "We now run the the solver and store the data at each time point to the initialized files. We also integrate A over the cytosolic volume at each time step to monitor the elevation in A over time and we display the concentration of A in the cytosol at t = 0.2 s."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "31d42278",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Set loglevel to warning in order not to pollute notebook output\n",
+ "logger.setLevel(logging.WARNING)\n",
+ "\n",
+ "avg_A = [A.initial_condition]\n",
+ "# define integration measure and total volume for computing average A at each time point\n",
+ "dx = d.Measure(\"dx\", domain=model_cur.cc['Cyto'].dolfin_mesh)\n",
+ "volume = d.assemble_mixed(1.0*dx)\n",
+ "# Solve\n",
+ "displayed = False\n",
+ "while model_cur.t < model_cur.final_t:\n",
+ " # Solve the system\n",
+ " model_cur.monolithic_solve()\n",
+ " # Update PSOCE\n",
+ " PSOCE.dolfin_function.vector()[:] = (d0/(dist_vals+d0)) * (1/(1+(AER.u[\"u\"].vector()[AER_idx]/cref)**4))\n",
+ " PSOCE.dolfin_function.vector().apply(\"insert\")\n",
+ " # Save results for post processing\n",
+ " for species_name, species in model_cur.sc.items:\n",
+ " results[species_name].write(model_cur.sc[species_name].u[\"u\"], model_cur.t)\n",
+ " SOCEfile.write(PSOCE.dolfin_function, model_cur.t)\n",
+ " # compute average A concentration at each time step\n",
+ " int_val = d.assemble_mixed(model_cur.sc['A'].u['u']*dx)\n",
+ " avg_A.append(int_val / volume)\n",
+ " if model_cur.t >= 0.2 and not displayed:\n",
+ " visualization.plot(model_cur.sc['A'].u['u'],\n",
+ " clip_plane=(1, 1, 0), clip_origin=(0.5, 0.5, 0.5))\n",
+ " displayed = True"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "621600d7",
+ "metadata": {},
+ "source": [
+ "Finally, we plot the average concentration of A over time."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b8651e71",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.plot(model_cur.tvec, avg_A, label='PM-ER gap = 50 nm')\n",
+ "other_case = np.loadtxt('/root/shared/gitrepos/smart-dev/examples/example5/results_ER_PM_gap_0.2/Avals.txt')\n",
+ "plt.plot(other_case[0,:], other_case[1,:], label='PM-ER gap = 200 nm')\n",
+ "plt.xlabel('Time (s)')\n",
+ "plt.ylabel('Cytosolic concentration of A (μM)')\n",
+ "plt.legend()\n",
+ "np.savetxt(str(result_folder / f\"Avals.txt\"), [model_cur.tvec, avg_A])"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "d302bf3f",
+ "metadata": {},
+ "source": [
+ "We also compare the AUC for A with previous numerical simulations (regression test)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "f8b3c866",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tvec = np.zeros(len(model_cur.tvec))\n",
+ "for i in range(len(model_cur.tvec)):\n",
+ " tvec[i] = float(model_cur.tvec[i])\n",
+ "auc_cur = np.trapz(np.array(avg_A), tvec)\n",
+ "auc_compare = 4.646230684534995\n",
+ "percent_error = 100*np.abs(auc_cur - auc_compare)/auc_compare\n",
+ "assert percent_error < .01,\\\n",
+ " f\"Failed regression test: Example 5 results deviate {percent_error:.3f}% from the previous numerical solution\""
+ ]
+ }
+ ],
+ "metadata": {
+ "jupytext": {
+ "cell_metadata_filter": "-all",
+ "main_language": "python",
+ "notebook_metadata_filter": "-all"
+ },
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.12"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/examples/example7/example7.ipynb b/examples/example7/example7.ipynb
new file mode 100644
index 0000000..2c86caa
--- /dev/null
+++ b/examples/example7/example7.ipynb
@@ -0,0 +1,304 @@
+{
+ "cells": [
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "f65f18d7",
+ "metadata": {},
+ "source": [
+ "# Example 7: Reaction-diffusion of molecule in a multicellular network\n",
+ "\n",
+ "Here, we consider the diffusion of a molecule from a central cell in a multicellular mesh to other cells in the environment."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "cc398816",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import dolfin as d\n",
+ "import sympy as sym\n",
+ "import numpy as np\n",
+ "import pathlib\n",
+ "import logging\n",
+ "import gmsh # must be imported before pyvista if dolfin is imported first\n",
+ "\n",
+ "from smart import config, mesh, model, mesh_tools, visualization\n",
+ "from smart.units import unit\n",
+ "from smart.model_assembly import (\n",
+ " Compartment,\n",
+ " Parameter,\n",
+ " Reaction,\n",
+ " Species,\n",
+ " SpeciesContainer,\n",
+ " ParameterContainer,\n",
+ " CompartmentContainer,\n",
+ " ReactionContainer,\n",
+ ")\n",
+ "\n",
+ "from matplotlib import pyplot as plt\n",
+ "import matplotlib.image as mpimg\n",
+ "from matplotlib import rcParams\n",
+ "\n",
+ "logger = logging.getLogger(\"smart\")\n",
+ "logger.setLevel(logging.INFO)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "95b9d865",
+ "metadata": {},
+ "source": [
+ "We define the relevant units here."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4f4023cf",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Aliases - base units\n",
+ "uM = unit.uM\n",
+ "um = unit.um\n",
+ "molecule = unit.molecule\n",
+ "sec = unit.sec\n",
+ "dimensionless = unit.dimensionless\n",
+ "# Aliases - units used in model\n",
+ "D_unit = um**2 / sec\n",
+ "flux_unit = uM * um / sec\n",
+ "vol_unit = uM\n",
+ "surf_unit = molecule / um**2"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "46582d26",
+ "metadata": {},
+ "source": [
+ "## Model generation\n",
+ "\n",
+ "We define the compartments and species first, with their respective containers."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "02a000f2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "EC = Compartment(\"EC\", 3, um, 1)\n",
+ "PM1 = Compartment(\"PM1\", 2, um, 11) # source cell\n",
+ "PM2 = Compartment(\"PM2\", 2, um, 12) # other cells\n",
+ "PM1.specify_nonadjacency(['PM2'])\n",
+ "PM2.specify_nonadjacency(['PM1'])\n",
+ "\n",
+ "cc = CompartmentContainer()\n",
+ "cc.add([EC, PM1, PM2])\n",
+ "\n",
+ "G = Species(\"G\", 1.0, vol_unit, 1000.0, D_unit, \"EC\")\n",
+ "Rbound = Species(\"Rbound\", 0.0, surf_unit, 0.0, D_unit, \"PM2\")\n",
+ "sc = SpeciesContainer()\n",
+ "sc.add([G, Rbound])"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "3c56e840",
+ "metadata": {},
+ "source": [
+ "Define parameters and reactions, then place in respective containers. Here, there are 3 reactions:\n",
+ "* r1: release of G from PM1\n",
+ "* r2: binding of G to PM2\n",
+ "* r3: degradation of G"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2e1f6882",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# G release from cell 1\n",
+ "j1pulse = Parameter(\"j1pulse\", 1000.0, flux_unit)\n",
+ "r1 = Reaction(\n",
+ " \"r1\",\n",
+ " [],\n",
+ " [\"G\"],\n",
+ " param_map={\"J\": \"j1pulse\"},\n",
+ " eqn_f_str=\"J\",\n",
+ " explicit_restriction_to_domain=\"PM1\",\n",
+ ")\n",
+ "# G binding to PM2\n",
+ "Rtot = Parameter(\"Rtot\", 100.0, surf_unit)\n",
+ "kbind = Parameter(\"kbind\", 1.0, 1/(uM*sec))\n",
+ "kunbind = Parameter(\"kunbind\", 0.01, 1/sec)\n",
+ "r2 = Reaction(\"r2\", [\"G\"], [\"Rbound\"],\n",
+ " param_map={\"on\":\"kbind\",\"off\":\"kunbind\",\"Rtot\":\"Rtot\"},\n",
+ " eqn_f_str=\"on*G*(Rtot-Rbound) - off*Rbound\",\n",
+ " explicit_restriction_to_domain=\"PM2\")\n",
+ "\n",
+ "# G degradation\n",
+ "kdeg = Parameter(\"kdeg\", 0.01, 1/sec)\n",
+ "r3 = Reaction(\"r3\", [\"G\"], [], param_map={\"k\":\"kdeg\"},\n",
+ " eqn_f_str=\"k*G\")\n",
+ "\n",
+ "pc = ParameterContainer()\n",
+ "pc.add([j1pulse,Rtot,kbind,kunbind,kdeg])\n",
+ "rc = ReactionContainer()\n",
+ "rc.add([r1,r2,r3])"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "15c35d39",
+ "metadata": {},
+ "source": [
+ "## Create and load in mesh\n",
+ "\n",
+ "Here, we consider cells embedded in a cube mesh. The source cell is located at (0,0,0) and 8 other cells are spread equidistant through the mesh."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "fe56e162",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "domain, facet_markers, cell_markers = mesh_tools.create_multicell(cubeSize=100.0, locVec1=[[0,0,0]],\n",
+ " locVec2=[[-30,-30,-30], [-30,30,-30], [-30,-30,30], [-30,30,30], \n",
+ " [30,-30,-30], [30,30,-30], [30,-30,30], [30,30,30]],\n",
+ " cellRad1 = 10.0, cellRad2 = 10.0, hCube = 5.0, hCell = 2.0)\n",
+ "# Write mesh and meshfunctions to file\n",
+ "mesh_folder = pathlib.Path(\"mesh\")\n",
+ "mesh_folder.mkdir(exist_ok=True)\n",
+ "mesh_path = mesh_folder / \"multicell_mesh.h5\"\n",
+ "mesh_tools.write_mesh(\n",
+ " domain, facet_markers, cell_markers, filename=mesh_path\n",
+ ")\n",
+ "parent_mesh = mesh.ParentMesh(\n",
+ " mesh_filename=str(mesh_path),\n",
+ " mesh_filetype=\"hdf5\",\n",
+ " name=\"parent_mesh\",\n",
+ ")\n",
+ "visualization.plot_dolfin_mesh(domain, cell_markers, facet_markers)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "0943588e",
+ "metadata": {},
+ "source": [
+ "Initialize model and solver."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ac88bdec",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "config_cur = config.Config()\n",
+ "config_cur.flags.update({\"allow_unused_components\": True})\n",
+ "model_cur = model.Model(pc, sc, cc, rc, config_cur, parent_mesh)\n",
+ "config_cur.solver.update(\n",
+ " {\n",
+ " \"final_t\": 10.0,\n",
+ " \"initial_dt\": 0.001,\n",
+ " \"time_precision\": 8,\n",
+ " \"reset_timestep_for_negative_solution\": True,\n",
+ " }\n",
+ ")\n",
+ "model_cur.initialize()"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "id": "5d5aacbd",
+ "metadata": {},
+ "source": [
+ "Initialize XDMF files for saving results, save model information to .pkl file, then solve the system until `model_cur.t > model_cur.final_t`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b54d28ca",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Write initial condition(s) to file\n",
+ "results = dict()\n",
+ "result_folder = pathlib.Path(f\"results\")\n",
+ "result_folder.mkdir(exist_ok=True)\n",
+ "for species_name, species in model_cur.sc.items:\n",
+ " results[species_name] = d.XDMFFile(\n",
+ " model_cur.mpi_comm_world, str(result_folder / f\"{species_name}.xdmf\")\n",
+ " )\n",
+ " results[species_name].parameters[\"flush_output\"] = True\n",
+ " results[species_name].write(model_cur.sc[species_name].u[\"u\"], model_cur.t)\n",
+ "model_cur.to_pickle(\"model_cur.pkl\")\n",
+ "\n",
+ "# Set loglevel to warning in order not to pollute notebook output\n",
+ "logger.setLevel(logging.WARNING)\n",
+ "# Solve\n",
+ "displayed = False\n",
+ "while True:\n",
+ " # Solve the system\n",
+ " model_cur.monolithic_solve()\n",
+ " model_cur.adjust_dt()\n",
+ " # Save results for post processing\n",
+ " for species_name, species in model_cur.sc.items:\n",
+ " results[species_name].write(model_cur.sc[species_name].u[\"u\"], model_cur.t)\n",
+ " print(f\"Done with t={model_cur.t}\")\n",
+ " # End if we've passed the final time\n",
+ " if model_cur.t >= model_cur.final_t:\n",
+ " break"
+ ]
+ }
+ ],
+ "metadata": {
+ "jupytext": {
+ "cell_metadata_filter": "-all",
+ "main_language": "python",
+ "notebook_metadata_filter": "-all"
+ },
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.12"
+ },
+ "vscode": {
+ "interpreter": {
+ "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
+ }
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/pyproject.toml b/pyproject.toml
index 0bbb4e3..e444256 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -3,7 +3,7 @@ requires = ["setuptools>=61.0.0", "wheel"]
[project]
name = "fenics-smart"
-version = "2.2.3"
+version = "2.3.0.beta.1"
description = "Spatial Modeling Algorithms for Reactions and Transport (SMART) is a high-performance finite-element-based simulation package for model specification and numerical simulation of spatially-varying reaction-transport processes in biological cells."
authors = [{ name = "Justin Laughlin", email = "justinglaughlin@gmail.com" }]
license = { file = "LICENSE" }
diff --git a/smart/mesh.py b/smart/mesh.py
index 96774ce..ba93abd 100644
--- a/smart/mesh.py
+++ b/smart/mesh.py
@@ -152,6 +152,10 @@ class ParentMesh(_Mesh):
as an xml or xdmf file containing a vertex mesh function of type double
extra_keys (list): list of names of extra keys to load from hdf5 mesh
file, specifying other subdomains besides main compartments (optional)
+ key_for_cells (string): a string specifying the extra key to use as a
+ mesh function to define cell markers in mesh (optional)
+ key_for_facets (string): a string specifying the extra key to use as a
+ mesh function to define facet markers in mesh (optional)
"""
mesh_filename: str
@@ -161,6 +165,8 @@ class ParentMesh(_Mesh):
use_partition: bool
curvature: d.MeshFunction
extra_keys: list = []
+ key_for_cells: str = ""
+ key_for_facets: str = ""
def __init__(
self,
@@ -171,6 +177,8 @@ def __init__(
mpi_comm=d.MPI.comm_world,
curvature=None,
extra_keys=[],
+ key_for_cells="",
+ key_for_facets="",
):
super().__init__(name)
self.use_partition = use_partition
@@ -202,6 +210,30 @@ def __init__(
# Otherwise just take what we got
self.curvature = curvature
+ # set cells and/or facets according to extra keys if applicable
+ if key_for_cells != "":
+ try:
+ idx = self.extra_keys.index(key_for_cells)
+ except ValueError:
+ raise ValueError(f"'{key_for_cells}' does not match an extra key")
+ assert self.subdomains[idx].dim() == self.dimensionality, (
+ f"Mesh function associated with '{key_for_cells}' "
+ "does not match mesh cell dimension"
+ )
+ self.mf["cells"] = self.subdomains[idx]
+ logger.info(f"Cell mesh function loaded from key '{key_for_cells}'")
+ if key_for_facets != "":
+ try:
+ idx = self.extra_keys.index(key_for_facets)
+ except ValueError:
+ raise ValueError(f"'{key_for_cells}' does not match an extra key")
+ assert self.subdomains[idx].dim() == self.dimensionality - 1, (
+ f"Mesh function associated with '{key_for_facets}' "
+ "does not match mesh facet dimension"
+ )
+ self.mf["facets"] = self.subdomains[idx]
+ logger.info(f"Facet mesh function loaded from key '{key_for_facets}'")
+
def get_mesh_from_id(self, id):
"Find the mesh that has the matching id."
# find the mesh in that has the matching id
@@ -290,8 +322,9 @@ def read_parent_mesh_functions_from_file(self):
assert len(self.child_meshes) > 0
# Init mesh functions
- self.mf["cells"] = self._read_parent_mesh_function_from_file(volume_dim)
- if self.has_surface:
+ if "cells" not in self.mf.keys():
+ self.mf["cells"] = self._read_parent_mesh_function_from_file(volume_dim)
+ if self.has_surface and "facets" not in self.mf.keys():
self.mf["facets"] = self._read_parent_mesh_function_from_file(surface_dim)
# If any cell markers are given as a list we also create mesh
diff --git a/smart/mesh_tools.py b/smart/mesh_tools.py
index e1a700a..4ffabef 100644
--- a/smart/mesh_tools.py
+++ b/smart/mesh_tools.py
@@ -719,6 +719,192 @@ def meshSizeCallback(dim, tag, x, y, z, lc):
return (dmesh, mf2, mf3)
+def create_multicell(
+ cubeSize: float = 100.0,
+ locVec1: list = [[0, 0, 0]],
+ locVec2: list = [],
+ cellRad1: float = 10.0,
+ cellRad2: float = 10.0,
+ hCube: float = 0,
+ hCell: float = 0,
+ interface_marker1: int = 11,
+ interface_marker2: int = 12,
+ outer_marker: int = 10,
+ extracell_tag: int = 1,
+ cell_vol_tag1: int = 2,
+ cell_vol_tag2: int = 3,
+ comm: MPI.Comm = d.MPI.comm_world,
+ verbose: bool = False,
+) -> Tuple[d.Mesh, d.MeshFunction, d.MeshFunction]:
+ """
+ Creates a mesh with an outer cube containing embedded cells at specified locations.
+ Args:
+ cubeSize: Length of cube sides
+ locVec: vector of cell locations
+ hCube: maximum mesh size for cube
+ hCell: maximum mesh size for cell surfaces
+ interface_marker: The value to mark facets on the interface with
+ outer_marker: The value to mark facets on the outer ellipsoid with
+ inner_vol_tag: The value to mark the inner spherical volume with
+ outer_vol_tag: The value to mark the outer spherical volume with
+ comm: MPI communicator to create the mesh with
+ verbose: If true print gmsh output, else skip
+ Returns:
+ A triplet (mesh, facet_marker, cell_marker)
+ """
+ import gmsh
+
+ if np.isclose(cubeSize, 0):
+ ValueError("Outer cube size is equal to zero")
+ if np.isclose(hCube, 0):
+ hCube = 0.1 * max(cubeSize)
+ if np.isclose(hCell, 0):
+ hCell = 0.2 * cubeSize if np.isclose(cellRad1, 0) else 0.2 * cellRad1
+ # if innerRad > outerRad or innerLength >= outerLength:
+ # ValueError("Inner cylinder does not fit inside outer cylinder")
+ # Create the two cylinder mesh using gmsh
+ gmsh.initialize()
+ gmsh.option.setNumber("General.Terminal", int(verbose))
+
+ gmsh.model.add("multicell")
+ # first add outer cube
+ cube = gmsh.model.occ.addBox(
+ -cubeSize / 2, -cubeSize / 2, -cubeSize / 2, cubeSize, cubeSize, cubeSize
+ )
+ if (np.isclose(cellRad1, 0) or len(locVec1) == 0) and (
+ np.isclose(cellRad2, 0) or len(locVec2) == 0
+ ):
+ # Just a cube!
+ gmsh.model.occ.synchronize()
+ gmsh.model.add_physical_group(3, [cube], tag=extracell_tag)
+ facets = gmsh.model.getBoundary([(3, cube)])
+ gmsh.model.add_physical_group(2, [facets[0][1]], tag=outer_marker)
+ else:
+ # Add cells
+ cell_list = []
+ # first add source cell(s)
+ for i in range(len(locVec1)):
+ cur_tag = gmsh.model.occ.addSphere(
+ locVec1[i][0], locVec1[i][1], locVec1[i][2], cellRad1
+ )
+ cell_list.append((3, cur_tag))
+ # now add additional cells
+ for i in range(len(locVec2)):
+ cur_tag = gmsh.model.occ.addSphere(
+ locVec2[i][0], locVec2[i][1], locVec2[i][2], cellRad2
+ )
+ cell_list.append((3, cur_tag))
+ # Create interface between cells and extracell
+ full_geo, maps = gmsh.model.occ.fragment([(3, cube)], cell_list)
+ cube_map = maps[0]
+ cell_maps = maps[1:]
+ gmsh.model.occ.synchronize()
+
+ # Get the outer boundary
+ outer_shells = gmsh.model.getBoundary(full_geo, oriented=False)
+ # Get the inner boundary
+ inner_shells1 = []
+ inner_shells2 = []
+ for i in range(0, len(locVec1)):
+ inner_shells1.append(gmsh.model.getBoundary(cell_maps[i], oriented=False))
+ for i in range(len(locVec1), len(cell_maps)):
+ inner_shells2.append(gmsh.model.getBoundary(cell_maps[i], oriented=False))
+ # Add physical markers for facets
+ gmsh.model.add_physical_group(2, [faces[1] for faces in outer_shells], tag=outer_marker)
+ gmsh.model.add_physical_group(
+ 2, [faces[0][1] for faces in inner_shells1], tag=interface_marker1
+ )
+ gmsh.model.add_physical_group(
+ 2, [faces[0][1] for faces in inner_shells2], tag=interface_marker2
+ )
+
+ # Physical markers for
+ all_volumes = [tag[1] for tag in cube_map]
+ inner_volumes1 = [tag[0][1] for tag in cell_maps[0 : len(locVec1)]]
+ inner_volumes2 = [tag[0][1] for tag in cell_maps[len(locVec1) :]]
+ outer_volume = []
+ for vol in all_volumes:
+ if (vol not in inner_volumes1) and (vol not in inner_volumes2):
+ outer_volume.append(vol)
+ gmsh.model.add_physical_group(3, outer_volume, tag=extracell_tag)
+ gmsh.model.add_physical_group(3, inner_volumes1, tag=cell_vol_tag1)
+ gmsh.model.add_physical_group(3, inner_volumes2, tag=cell_vol_tag2)
+
+ def meshSizeCallback(dim, tag, x, y, z, lc):
+ # mesh length is hEdge at the PM (defaults to 0.1*outerRad,
+ # or set when calling function) and hInnerEdge at the ERM
+ # (defaults to 0.2*innerRad, or set when calling function)
+ # between these, the value is interpolated based on r (polar coord),
+ # and inside the value is interpolated between hInnerEdge and 0.2*innerRad
+ # if innerRad=0, then the mesh length is interpolated between
+ # hEdge at the PM and 0.2*outerRad in the center
+
+ if (np.isclose(cellRad1, 0) or len(locVec1) == 0) and (
+ np.isclose(cellRad2, 0) or len(locVec2) == 0
+ ):
+ return hCube
+ elif np.isclose(cellRad1, 0) or len(locVec1) == 0:
+ cell_locs1 = [np.inf]
+ cell_locs2 = np.sqrt(
+ (x - np.array(locVec2)[:, 0]) ** 2
+ + (y - np.array(locVec2)[:, 1]) ** 2
+ + (z - np.array(locVec2)[:, 2]) ** 2
+ )
+ elif np.isclose(cellRad2, 0) or len(locVec2) == 0:
+ cell_locs1 = np.sqrt(
+ (x - np.array(locVec1)[:, 0]) ** 2
+ + (y - np.array(locVec1)[:, 1]) ** 2
+ + (z - np.array(locVec1)[:, 2]) ** 2
+ )
+ cell_locs2 = [np.inf]
+ else:
+ cell_locs1 = np.sqrt(
+ (x - np.array(locVec1)[:, 0]) ** 2
+ + (y - np.array(locVec1)[:, 1]) ** 2
+ + (z - np.array(locVec1)[:, 2]) ** 2
+ )
+ cell_locs2 = np.sqrt(
+ (x - np.array(locVec2)[:, 0]) ** 2
+ + (y - np.array(locVec2)[:, 1]) ** 2
+ + (z - np.array(locVec2)[:, 2]) ** 2
+ )
+ closest_cell1 = min(cell_locs1)
+ closest_cell2 = min(cell_locs2)
+ if (closest_cell1 < cellRad1) or (closest_cell2 < cellRad2):
+ return hCell
+ else:
+ if closest_cell1 < closest_cell2:
+ cellWeight = np.exp(-(closest_cell1 - cellRad1) / (0.2 * cellRad1))
+ else:
+ cellWeight = np.exp(-(closest_cell2 - cellRad2) / (0.2 * cellRad2))
+ return hCell * cellWeight + hCube * (1 - cellWeight)
+
+ gmsh.model.mesh.setSizeCallback(meshSizeCallback)
+ # set off the other options for mesh size determination
+ gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
+ gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
+ gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)
+ # this changes the algorithm from Frontal-Delaunay to Delaunay,
+ # which may provide better results when there are larger gradients in mesh size
+ gmsh.option.setNumber("Mesh.Algorithm", 5)
+
+ gmsh.model.mesh.generate(3)
+ rank = MPI.COMM_WORLD.rank
+ tmp_folder = pathlib.Path(f"tmp_extracell_{cubeSize}_{cellRad1}_{cellRad2}_{rank}")
+ tmp_folder.mkdir(exist_ok=True)
+ gmsh_file = tmp_folder / "extracell.msh"
+ gmsh.write(str(gmsh_file))
+ gmsh.finalize()
+
+ # return dolfin mesh of max dimension (parent mesh) and marker functions mf2 and mf3
+ dmesh, mf2, mf3 = gmsh_to_dolfin(str(gmsh_file), tmp_folder, 3, comm)
+ # remove tmp mesh and tmp folder
+ gmsh_file.unlink(missing_ok=False)
+ tmp_folder.rmdir()
+ # return dolfin mesh, mf2 (2d tags) and mf3 (3d tags)
+ return (dmesh, mf2, mf3)
+
+
def create_ellipses(
xrad_outer: float = 3.0,
yrad_outer: float = 1.0,
@@ -1300,6 +1486,291 @@ def meshSizeCallback(dim, tag, x, y, z, lc):
return (dmesh, mf2, mf3)
+def create_2Dcell_xy(
+ outerExpr: str = "",
+ innerExpr: str = "",
+ hEdge: float = 0,
+ hInnerEdge: float = 0,
+ interface_marker: int = 12,
+ outer_marker: int = 10,
+ inner_tag: int = 2,
+ outer_tag: int = 1,
+ comm: MPI.Comm = d.MPI.comm_world,
+ verbose: bool = False,
+ half_cell: bool = True,
+ return_curvature: bool = False,
+) -> Tuple[d.Mesh, d.MeshFunction, d.MeshFunction]:
+ """
+ Creates a 2D mesh of a cell profile, with the bounding curve defined in
+ terms of r and z (e.g. unit circle would be "r**2 + (z-1)**2 - 1)
+ It is assumed that substrate is present at z = 0, so if the curve extends
+ below z = 0 , there is a sharp cutoff.
+ If half_cell = True, only have of the contour is constructed, with a
+ left zero-flux boundary at r = 0.
+ Can include one compartment inside another compartment.
+ Recommended for use with the axisymmetric feature of SMART.
+
+ Args:
+ outerExpr: String implicitly defining an r-z curve for the outer surface
+ innerExpr: String implicitly defining an r-z curve for the inner surface
+ hEdge: maximum mesh size at the outer edge
+ hInnerEdge: maximum mesh size at the edge
+ of the inner compartment
+ interface_marker: The value to mark facets on the interface with
+ outer_marker: The value to mark facets on edge of the outer ellipse with
+ inner_tag: The value to mark the inner ellipse surface with
+ outer_tag: The value to mark the outer ellipse surface with
+ comm: MPI communicator to create the mesh with
+ verbose: If true print gmsh output, else skip
+ half_cell: If true, consider r=0 the symmetry axis for an axisymm shape
+ Returns:
+ Tuple (mesh, facet_marker, cell_marker)
+ """
+ import gmsh
+
+ if outerExpr == "":
+ ValueError("Outer surface is not defined")
+
+ if return_curvature:
+ # create full mesh for curvature analysis and then map onto half mesh
+ # if half_cell_with_curvature is True
+ half_cell_with_curvature = half_cell
+ half_cell = False
+
+ rValsOuter, zValsOuter = implicit_curve(outerExpr)
+
+ if not innerExpr == "":
+ rValsInner, zValsInner = implicit_curve(innerExpr)
+ zMid = np.mean(zValsInner)
+ ROuterVec = np.sqrt(rValsOuter**2 + (zValsOuter - zMid) ** 2)
+ RInnerVec = np.sqrt(rValsInner**2 + (zValsInner - zMid) ** 2)
+ maxOuterDim = max(ROuterVec)
+ maxInnerDim = max(RInnerVec)
+ else:
+ zMid = np.mean(zValsOuter)
+ ROuterVec = np.sqrt(rValsOuter**2 + (zValsOuter - zMid) ** 2)
+ maxOuterDim = max(ROuterVec)
+ if np.isclose(hEdge, 0):
+ hEdge = 0.1 * maxOuterDim
+ if np.isclose(hInnerEdge, 0):
+ hInnerEdge = 0.2 * maxOuterDim if innerExpr == "" else 0.2 * maxInnerDim
+ # Create the 2D mesh using gmsh
+ gmsh.initialize()
+ gmsh.option.setNumber("General.Terminal", int(verbose))
+ gmsh.model.add("2DCell")
+ # first add outer body and revolve
+ outer_tag_list = []
+ for i in range(len(rValsOuter)):
+ cur_tag = gmsh.model.occ.add_point(rValsOuter[i], zValsOuter[i], 0.0)
+ outer_tag_list.append(cur_tag)
+ outer_spline_tag = gmsh.model.occ.add_spline(outer_tag_list)
+ if not half_cell:
+ outer_tag_list2 = []
+ for i in range(len(rValsOuter)):
+ cur_tag = gmsh.model.occ.add_point(-rValsOuter[i], zValsOuter[i], 0.0)
+ outer_tag_list2.append(cur_tag)
+ outer_spline_tag2 = gmsh.model.occ.add_spline(outer_tag_list2)
+ if np.isclose(zValsOuter[-1], 0): # then include substrate at z=0
+ if half_cell:
+ origin_tag = gmsh.model.occ.add_point(0, 0, 0)
+ symm_axis_tag = gmsh.model.occ.add_line(origin_tag, outer_tag_list[0])
+ bottom_tag = gmsh.model.occ.add_line(origin_tag, outer_tag_list[-1])
+ outer_loop_tag = gmsh.model.occ.add_curve_loop(
+ [outer_spline_tag, bottom_tag, symm_axis_tag]
+ )
+ else:
+ bottom_tag = gmsh.model.occ.add_line(outer_tag_list[-1], outer_tag_list2[-1])
+ outer_loop_tag = gmsh.model.occ.add_curve_loop(
+ [outer_spline_tag, outer_spline_tag2, bottom_tag]
+ )
+ else:
+ if half_cell:
+ symm_axis_tag = gmsh.model.occ.add_line(outer_tag_list[0], outer_tag_list[-1])
+ outer_loop_tag = gmsh.model.occ.add_curve_loop([outer_spline_tag, symm_axis_tag])
+ else:
+ outer_loop_tag = gmsh.model.occ.add_curve_loop([outer_spline_tag, outer_spline_tag2])
+ cell_plane_tag = gmsh.model.occ.add_plane_surface([outer_loop_tag])
+
+ if innerExpr == "":
+ # No inner shape in this case
+ gmsh.model.occ.synchronize()
+ gmsh.model.add_physical_group(2, [cell_plane_tag], tag=outer_tag)
+ facets = gmsh.model.getBoundary([(2, cell_plane_tag)])
+ facet_tag_list = []
+ for i in range(len(facets)):
+ facet_tag_list.append(facets[i][1])
+ if half_cell: # if half, set symmetry axis to 0 (no flux)
+ xmin, ymin, zmin = (-hInnerEdge / 10, -1, -hInnerEdge / 10)
+ xmax, ymax, zmax = (hInnerEdge / 10, max(zValsOuter) + 1, hInnerEdge / 10)
+ all_symm_bound = gmsh.model.occ.get_entities_in_bounding_box(
+ xmin, ymin, zmin, xmax, ymax, zmax, dim=1
+ )
+ symm_bound_markers = []
+ for i in range(len(all_symm_bound)):
+ symm_bound_markers.append(all_symm_bound[i][1])
+ gmsh.model.add_physical_group(1, symm_bound_markers, tag=0)
+ gmsh.model.add_physical_group(1, facet_tag_list, tag=outer_marker)
+ else:
+ # Add inner shape
+ inner_tag_list = []
+ for i in range(len(rValsInner)):
+ cur_tag = gmsh.model.occ.add_point(rValsInner[i], zValsInner[i], 0.0)
+ inner_tag_list.append(cur_tag)
+ inner_spline_tag = gmsh.model.occ.add_spline(inner_tag_list)
+ if half_cell:
+ symm_inner_tag = gmsh.model.occ.add_line(inner_tag_list[0], inner_tag_list[-1])
+ inner_loop_tag = gmsh.model.occ.add_curve_loop([inner_spline_tag, symm_inner_tag])
+ else:
+ inner_tag_list2 = []
+ for i in range(len(rValsInner)):
+ cur_tag = gmsh.model.occ.add_point(-rValsInner[i], zValsInner[i], 0.0)
+ inner_tag_list2.append(cur_tag)
+ inner_spline_tag2 = gmsh.model.occ.add_spline(inner_tag_list2)
+ inner_loop_tag = gmsh.model.occ.add_curve_loop([inner_spline_tag, inner_spline_tag2])
+ inner_plane_tag = gmsh.model.occ.add_plane_surface([inner_loop_tag])
+ cell_plane_list = [cell_plane_tag]
+ inner_plane_list = [inner_plane_tag]
+
+ outer_volume = []
+ inner_volume = []
+ all_volumes = []
+ inner_marker_list = []
+ outer_marker_list = []
+ for i in range(len(cell_plane_list)):
+ cell_plane_tag = cell_plane_list[i]
+ inner_plane_tag = inner_plane_list[i]
+ # Create interface between 2 objects
+ two_shapes, (outer_shape_map, inner_shape_map) = gmsh.model.occ.fragment(
+ [(2, cell_plane_tag)], [(2, inner_plane_tag)]
+ )
+ gmsh.model.occ.synchronize()
+
+ # Get the outer boundary
+ outer_shell = gmsh.model.getBoundary(two_shapes, oriented=False)
+ for i in range(len(outer_shell)):
+ outer_marker_list.append(outer_shell[i][1])
+ # Get the inner boundary
+ inner_shell = gmsh.model.getBoundary(inner_shape_map, oriented=False)
+ for i in range(len(inner_shell)):
+ inner_marker_list.append(inner_shell[i][1])
+ for tag in outer_shape_map:
+ all_volumes.append(tag[1])
+ for tag in inner_shape_map:
+ inner_volume.append(tag[1])
+
+ for vol in all_volumes:
+ if vol not in inner_volume:
+ outer_volume.append(vol)
+
+ # Add physical markers for facets
+ if half_cell: # if half, set symmetry axis to 0 (no flux)
+ xmin, ymin, zmin = (-hInnerEdge / 10, -1, -hInnerEdge / 10)
+ xmax, ymax, zmax = (hInnerEdge / 10, max(zValsOuter) + 1, hInnerEdge / 10)
+ all_symm_bound = gmsh.model.occ.get_entities_in_bounding_box(
+ xmin, ymin, zmin, xmax, ymax, zmax, dim=1
+ )
+ symm_bound_markers = []
+ for i in range(len(all_symm_bound)):
+ symm_bound_markers.append(all_symm_bound[i][1])
+ # note that this first call sets the symmetry axis to tag 0 and
+ # this is not overwritten by the next calls to add_physical_group
+ gmsh.model.add_physical_group(1, symm_bound_markers, tag=0)
+ gmsh.model.add_physical_group(1, outer_marker_list, tag=outer_marker)
+ gmsh.model.add_physical_group(1, inner_marker_list, tag=interface_marker)
+
+ # Physical markers for "volumes"
+ gmsh.model.add_physical_group(2, outer_volume, tag=outer_tag)
+ gmsh.model.add_physical_group(2, inner_volume, tag=inner_tag)
+
+ def meshSizeCallback(dim, tag, x, y, z, lc):
+ # mesh length is hEdge at the PM and hInnerEdge at the inner membrane
+ # between these, the value is interpolated based on the relative distance
+ # between the two membranes.
+ # Inside the inner shape, the value is interpolated between hInnerEdge
+ # and lc3, where lc3 = max(hInnerEdge, 0.2*maxInnerDim)
+ # if innerRad=0, then the mesh length is interpolated between
+ # hEdge at the PM and 0.2*maxOuterDim in the center
+ rCur = np.sqrt(x**2 + z**2)
+ RCur = np.sqrt(rCur**2 + (y - zMid) ** 2)
+ outer_dist = np.sqrt((rCur - rValsOuter) ** 2 + (z - zValsOuter) ** 2)
+ np.append(outer_dist, z) # include the distance from the substrate
+ dist_to_outer = min(outer_dist)
+ if innerExpr == "":
+ lc3 = 0.2 * maxOuterDim
+ dist_to_inner = RCur
+ in_outer = True
+ else:
+ inner_dist = np.sqrt((rCur - rValsInner) ** 2 + (y - zValsInner) ** 2)
+ dist_to_inner = min(inner_dist)
+ inner_idx = np.argmin(inner_dist)
+ inner_rad = RInnerVec[inner_idx]
+ R_rel_inner = RCur / inner_rad
+ lc3 = max(hInnerEdge, 0.2 * maxInnerDim)
+ in_outer = R_rel_inner > 1
+ lc1 = hEdge
+ lc2 = hInnerEdge
+ if in_outer:
+ lcTest = lc1 + (lc2 - lc1) * (dist_to_outer) / (dist_to_inner + dist_to_outer)
+ else:
+ lcTest = lc2 + (lc3 - lc2) * (1 - R_rel_inner)
+ return lcTest
+
+ gmsh.model.mesh.setSizeCallback(meshSizeCallback)
+ # set off the other options for mesh size determination
+ gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
+ gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
+ gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)
+ # this changes the algorithm from Frontal-Delaunay to Delaunay,
+ # which may provide better results when there are larger gradients in mesh size
+ gmsh.option.setNumber("Mesh.Algorithm", 5)
+
+ gmsh.model.mesh.generate(2)
+ rank = MPI.COMM_WORLD.rank
+ tmp_folder = pathlib.Path(f"tmp_2DCell_{rank}")
+ tmp_folder.mkdir(exist_ok=True)
+ gmsh_file = tmp_folder / "2DCell.msh"
+ gmsh.write(str(gmsh_file))
+ gmsh.finalize()
+
+ # return dolfin mesh of max dimension (parent mesh) and marker functions mf2 and mf3
+ dmesh, mf1, mf2 = gmsh_to_dolfin(str(gmsh_file), tmp_folder, 2, comm)
+ # remove tmp mesh and tmp folder
+ gmsh_file.unlink(missing_ok=False)
+ tmp_folder.rmdir()
+ # return dolfin mesh, mf1 (1d tags) and mf2 (2d tags)
+ if return_curvature:
+ if innerExpr == "":
+ facet_list = [outer_marker]
+ cell_list = [outer_tag]
+ else:
+ facet_list = [outer_marker, interface_marker]
+ cell_list = [outer_tag, inner_tag]
+ if half_cell_with_curvature: # will likely not work in parallel...
+ dmesh_half, mf1_half, mf2_half = create_2Dcell(
+ outerExpr,
+ innerExpr,
+ hEdge,
+ hInnerEdge,
+ interface_marker,
+ outer_marker,
+ inner_tag,
+ outer_tag,
+ comm,
+ verbose,
+ half_cell=True,
+ return_curvature=False,
+ )
+ kappa_mf = compute_curvature(
+ dmesh, mf1, mf2, facet_list, cell_list, half_mesh_data=(dmesh_half, mf1_half)
+ )
+ (dmesh, mf1, mf2) = (dmesh_half, mf1_half, mf2_half)
+ else:
+ kappa_mf = compute_curvature(dmesh, mf1, mf2, facet_list, cell_list)
+ return (dmesh, mf1, mf2, kappa_mf)
+ else:
+ return (dmesh, mf1, mf2)
+
+
def gmsh_to_dolfin(
gmsh_file_name: str,
tmp_folder: pathlib.Path = pathlib.Path("tmp_folder"),
@@ -1340,6 +1811,8 @@ def gmsh_to_dolfin(
# convert cell mesh
cells = mesh_in.get_cells_type(cell_type)
cell_data = mesh_in.get_cell_data("gmsh:physical", cell_type) # extract values of tags
+ if dimension == 2 and np.all(mesh_in.points[:, 2] == 0):
+ mesh_in.points = mesh_in.points[:, :2] # then prune z values
out_mesh_cell = meshio.Mesh(
points=mesh_in.points,
cells={cell_type: cells},
diff --git a/smart/model.py b/smart/model.py
index 8f51475..737f911 100644
--- a/smart/model.py
+++ b/smart/model.py
@@ -1035,6 +1035,16 @@ def _init_4_7_set_initial_conditions(self):
new_vals = vec_new[mesh_map] # reorder to match dof ordering
vec.set_local(new_vals)
vec.apply("insert")
+ elif parameter.type == ParameterType.mesh_quantity:
+ if parameter.compartment not in self.cc.keys:
+ raise ValueError(
+ f"Compartment name {parameter.compartment} for parameter"
+ f"{parameter.name} does not match a known compartment"
+ )
+ V_cur = d.FunctionSpace(self.cc[parameter.compartment].dolfin_mesh, "P", 1)
+ parameter.dolfin_function = d.Function(V_cur)
+ parameter.dolfin_function.vector()[:] = parameter.value
+ parameter.dolfin_function.vector().apply("insert")
def _init_5_1_reactions_to_fluxes(self):
"""Convert reactions to flux objects"""
diff --git a/smart/model_assembly.py b/smart/model_assembly.py
index b8c6ef0..f3094ae 100644
--- a/smart/model_assembly.py
+++ b/smart/model_assembly.py
@@ -64,6 +64,7 @@ class ParameterType(str, Enum):
constant = "constant"
expression = "expression"
from_xdmf = "from_xdmf"
+ mesh_quantity = "mesh_quantity"
class InvalidObjectException(Exception):
@@ -602,6 +603,7 @@ class Parameter(ObjectInstance):
is_time_dependent: bool = False
is_space_dependent: bool = False
compartment: str = ""
+ mesh_quantity: bool = False
def to_dict(self):
"""Convert to a dict that can be used to recreate the object."""
@@ -731,6 +733,36 @@ def from_xdmf(
return parameter
+ @classmethod
+ def mesh_quantity(
+ cls, name, init_val, unit, compartment, group="", notes="", use_preintegration=False
+ ):
+ """ "
+ Initialize as a generic dolfin function over the mesh.
+ """
+ logger.debug(f"Initializing parameter {name} as mesh quantity")
+ if use_preintegration:
+ logger.warning(
+ f"Setting use_preintegration to False for parameter {name}."
+ "Not currently implemented for parameters given as mesh quantities"
+ )
+ use_preintegration = False
+ parameter = cls(
+ name,
+ init_val,
+ unit,
+ group=group,
+ notes=notes,
+ use_preintegration=use_preintegration,
+ )
+ parameter.compartment = compartment
+ # initialize instance
+ parameter.is_time_dependent = False
+ parameter.is_space_dependent = True
+ parameter.type = ParameterType.mesh_quantity
+ parameter.__post_init__()
+ return parameter
+
@classmethod
def from_expression(
cls,
@@ -846,7 +878,7 @@ def __post_init__(self):
@property
def dolfin_quantity(self):
- if self.type == ParameterType.from_xdmf:
+ if self.type == ParameterType.from_xdmf or self.type == ParameterType.mesh_quantity:
return self.dolfin_function * self.unit
elif hasattr(self, "dolfin_expression"):
return self.dolfin_expression * self.unit