Update for skin. Working version for binding.

content/skin/SkinDiffusion.ipynb

@@ -27,7 +27,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -38,7 +38,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -76,11 +76,11 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Loading Domain {gridName}...\n",
+      "Loading Domain 'skin2d-aniso.ugx'...\n",
       "Domain loaded.\n",
       "Refining ...\n",
-      "Refining step {i} ...\n",
-      "Refining step {i} ...\n",
+      "Refining step {0} ...\n",
+      "Refining step {1} ...\n",
       "Refining done\n"
      ]
     }

content/skin/SkinDiffusion.py

@@ -1,146 +1,260 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# [<img src="../../header.svg">](../index.ipynb)
+# 
+# ---
+
+# # Application: Simulation of Drug Transport across a Virtual Skin Membrane
+# 
+# [Advanced version](SkinDiffusionWithLagtime.ipybnd)
+
+# ### Preliminaries
+
+# In[11]:
+
+
 import ug4py.pyugcore as ug4
 import ug4py.pylimex as limex
 import ug4py.pyconvectiondiffusion as cd
-# import ug4py.pysuperlu as slu
 
 
-# Setup:
-# defining needed subsets, gird and number of refinements
+# In[12]:
+
+
+# Importing
+import sys
+sys.path.append('..')
+import modsimtools as util
+
+
+# ## Solve problem
+
+# In[3]:
+
+
+# Setup: Defining requires subsets, grid and number of refinements
 requiredSubsets = ["LIP", "COR", "BOTTOM_SC", "TOP_SC"] # defining subsets
-gridName = "skin2d-aniso.ugx"  # Grid created with ProMesh
+gridName = 'skin2d-aniso.ugx'  # Grid created with ProMesh
 numRefs = 2  # Number of Refinement steps on given grid
 
-# Choosing a domain object
-# (either 1d, 2d or 3d suffix)
-dom = ug4.Domain2d()
-
-# Loading the given grid into the domain
-print(f"Loading Domain {gridName} ...")
-ug4.LoadDomain(dom, gridName)
-print("Domain loaded.")
-
-# Optional: Refining the grid
-if numRefs > 0:
-    print("Refining ...")
-    refiner = ug4.GlobalDomainRefiner(dom)
-    for i in range(numRefs):
-        ug4.TerminateAbortedRun()
-        refiner.refine()
-        print(f"Refining step {i} ...")
-
-    print("Refining done")
-
-# checking if geometry has the needed subsets of the probelm
-sh = dom.subset_handler()
-for e in requiredSubsets:
-    if sh.get_subset_index(e) == -1:
-        print(f"Domain does not contain subset {e}.")
-        sys.exit(1)
+
+# In[4]:
+
+
+dom = util.CreateDomain(gridName, numRefs, requiredSubsets)
+
+
+# ### Create Approximation space
+
+# In[5]:
+
 
 # Create approximation space which describes the unknowns in the equation
-fct = "u"  # name of the function
-type = "Lagrange"
-order = 1  # polynom order for lagrange
-
-approxSpace = ug4.ApproximationSpace2d(dom)
-approxSpace.add_fct(fct, type, order)
-approxSpace.init_levels()
-approxSpace.init_surfaces()  
-approxSpace.init_top_surface()
-approxSpace.print_statistic()
-
-# Create discretization for the subsets LIP and COR
-# perform the discretization on the actual elements
-# using first order finite volumes (FV1) on the 2d grid
-KDcor=1.0
-KCor=1.0
+approxSpaceDesc = dict(fct = "u", type = "Lagrange", order = 1)
+
+
+# In[6]:
+
+
+approxSpace = util.CreateApproximationSpace(dom, approxSpaceDesc)
+
+
+# ## Create a convection-diffusion-equation
+# $$\frac{\partial Ku}{\partial t} + \nabla \cdot [-DK \nabla u] = 0$$
+# Partition coefficients
+
+# In[7]:
+
+
+K={}
+K["COR"]=1e-3
+K["LIP"]=1.0
+
+
+# Diffusion coefficients
+
+# In[8]:
+
+
+D={}
+D["COR"]=1.0
+D["LIP"]=1.0
+
+
+# Creating two instances of a convection diffusion equation
+
+# In[9]:
+
 
 elemDisc={}
-# creating instance of a convection diffusion equation
-elemDisc["COR"] = cd.ConvectionDiffusionFV12d("u", "COR")
-elemDisc["COR"].set_diffusion(KDcor)
-elemDisc["COR"].set_mass_scale(KCor)
+elemDisc["COR"] = util.CreateDiffusionElemDisc("u", "COR", K["COR"], K["COR"]*D["COR"], 0.0)
+elemDisc["LIP"] = util.CreateDiffusionElemDisc("u", "LIP", K["LIP"], K["LIP"]*D["LIP"], 0.0)
 
-elemDisc["LIP"] = cd.ConvectionDiffusionFV12d("u", "LIP")
-elemDisc["LIP"].set_diffusion(KDcor)
-elemDisc["LIP"].set_mass_scale(KCor)
 
+# ### Boundary conditions
 # ug4 separates the boundary value and the discretization
 # boundary conditions can be enforced through a post-process (dirichlet).
 # To init at boundary, the value, function name from the Approximationspace
-# and the subset name are needed
+# #and the subset name are needed
+
+# In[10]:
+
+
 dirichletBND = ug4.DirichletBoundary2dCPU1()
 dirichletBND.add(1.0, "u", "TOP_SC")
 dirichletBND.add(0.0, "u", "BOTTOM_SC")
 
 
-# create the discretization object which combines all the
-# separate discretizations into one domain discretization.
+# ### Global discretization
+# Create the discretization object which combines all the separate discretizations into one domain discretization.
+
+# In[11]:
+
+
 domainDisc = ug4.DomainDiscretization2dCPU1(approxSpace)
 domainDisc.add(elemDisc["COR"])
 domainDisc.add(elemDisc["LIP"])
 domainDisc.add(dirichletBND)
 
 
-# Create Solver
-# In this case we use an LU (Lower Upper) Solver for an exact solution
+# ## Create  solver
+
+# In[12]:
+
+
 lsolver=ug4.LUCPU1()
+
+
+# In[13]:
+
+
+# import pysuperlu as slu
 # lsolver=slu.SuperLUCPU1()
 
-# Solve the transient problem
-# Use the Approximationspace to
-# create a vector of unknowns and a vector
-# which contains the right hand side
-usol = ug4.GridFunction2dCPU1(approxSpace)
 
-# Init the vector representing the unknowns with 0 to obtain
-# reproducable results
-ug4.Interpolate(0.0, usol, "u")
+# ## Solve transient problem
+# 
+# Solve the transient problem.
+# Use the Approximationspace to create a vector of unknowns and a vector which contains the right hand side
+
+# In[14]:
+
 
 # Define start time, end time and step size
 startTime = 0.0
 endTime = 1000.0
 dt = 25.0
 
-# create a time discretization with the theta-scheme
-# takes in domain and a theta
-# theta = 1.0 -> Implicit Euler, 
-# 0.0 -> Explicit Euler 
+
+# Create a time discretization (with the theta-scheme).
+# Requires domain and theta ($\theta$ = 1.0 -> implicit Euler, 
+# $\theta$ = 0.0 -> explicit Euler )
+# 
+
+# In[15]:
+
+
 timeDisc=ug4.ThetaTimeStepCPU1(domainDisc, 1.0) 
 
-# creating Time Iterator, setting the solver and step size
-timeInt = limex.LinearTimeIntegrator2dCPU1(timeDisc)
-#timeInt = limex.ConstStepLinearTimeIntegrator2dCPU1(timeDisc)
+
+# Create time integrator (requires solver and step size).
+
+# In[16]:
+
+
+#timeInt = limex.LinearTimeIntegrator2dCPU1(timeDisc)
+timeInt = limex.ConstStepLinearTimeIntegrator2dCPU1(timeDisc)
 timeInt.set_linear_solver(lsolver)
 timeInt.set_time_step(dt)
 
 
+# Solving the transient problem
+
+# In[17]:
+
+
+# Create grid function for solution.
+usol = ug4.GridFunction2dCPU1(approxSpace)
+
+# Init the vector representing the unknowns with 0 to obtain
+# reproducable results
+ug4.Interpolate(0.0, usol, "u")
+
 # Exporting the result to a vtu-file
 # can be visualized in paraview or with a python extension
-# We add a corresponding observer to the time integrator.
-def MyCallback(usol, step, time, dt) :
-    ug4.WriteGridFunctionToVTK(usol, "vtk/SkinDiffusion_"+str(int(step)).zfill(5)+".vtu")
-pyobserver =ug4.PythonCallbackObserver2dCPU1(MyCallback) 
-timeInt.attach_observer(pyobserver)
+ug4.WriteGridFunctionToVTK(usol, "vtk/SkinDiffusion_Initial")
+
+
+# In[18]:
+
 
-# Solve the transient problem.
 try:
     timeInt.apply(usol, endTime, usol, startTime)
 except Exception as inst:
     print(inst)
 
 
+# In[19]:
+
+
+# Write last solution.
+ug4.WriteGridFunctionToVTK(usol, "vtk/SkinDiffusion_Final")
 
-# Plot the result (if Pyvista is present)
-hasPyvista = True
+
+# ## Plotting the result using pyvista
+
+# In[20]:
+
+
+hasPyvista=True
 try:
     import pyvista
 except Exception as inst:
-    hasPyvista = False
+    print(inst)
+    hasPyvista=False
+
+
+# In[21]:
+
+
+if hasPyvista:
+    pyvista.start_xvfb()
+    result = pyvista.read('vtk/SkinDiffusion_Final.vtu')
+
+    pyvista.set_jupyter_backend('static')
+    #pyvista.set_jupyter_backend('trame')
 
-if hasPyvista :
-    result = pyvista.read('vtk/SkinDiffusion_00040.vtu')
     print("Pyvista input: ")
-    print(result)
-    result.plot(scalars="u", show_edges=True, cmap='coolwarm')
-#scalars="node_value", categories=True
+    result.plot(scalars="u", show_edges=True, cmap='jet')
+
+
+# In[ ]:
+
+
+
+
+
+# In[ ]:
+
+
+
+
+
+# In[ ]:
+
+
+
+
+
+# In[ ]:
+
+
+
+
+
+# In[ ]:
+
+
+
+