This notebook shows how to describe a simple subset of $R^d$, and to compute finite differences taking into account a given Dirichlet data on the boundary. It is a prerequisite to the implementation of monotone schemes on non-rectangular domains.
Note: The library presented in this notebook features some elementary primitives of algorithmic geometry. However, their expressivity is rather limited, and exact predicates are not implemented. Therefore, please consider using a dedicated computational geometry package if you intend to go much beyond the basic examples presented in this series of notebooks.
Summary of volume Algorithmic tools, this series of notebooks.
Main summary of the Adaptive Grid Discretizations book of notebooks, including the other volumes.
Acknowledgement. Some of the experiments presented in these notebooks are part of ongoing research with Ludovic Métivier and Da Chen.
Copyright Jean-Marie Mirebeau, Centre Borelli, ENS Paris-Saclay, CNRS, University Paris-Saclay
import sys; sys.path.insert(0,"..") # Allow import of agd from parent directory (useless if conda package installed)
#from Miscellaneous import TocTools; print(TocTools.displayTOC('SubsetRd','Algo'))
from agd import Domain
from agd import AutomaticDifferentiation as ad
import scipy.linalg
import numpy as np
from matplotlib import pyplot as plt
def reload_packages():
import importlib
Domain = importlib.reload(sys.modules['agd.Domain'])
ad = importlib.reload(sys.modules['agd.AutomaticDifferentiation'])
ad.reload_submodules()
The provided library allows to define a family of basic shapes, combine them with boolean operations, and affine transformations. Once the desired shape is constructed, two main functions are available: a level set map, and the distance to the boundary along a given half line.
Define a ball and a box domain, here in dimension two, but these primitives are dimension independent.
ball = Domain.Ball([0,0]) # Two dimensional ball centered at the origin with radius one.
box = Domain.Box([[0,1],[0,1]]) # Square domain [0,1]x[0,1]
The standard boolean operations allow to combine elementary shapes: intersection, union, and relative complement.
cap = Domain.Intersection(ball,box)
cup = Domain.Union(ball,box)
compl = Domain.Complement(ball,box)
Bands and Convex polygons are some potentially useful additional shapes:
band = Domain.Band([1.,2.],[-2.,1.]) # Direction, bounds
triangle = Domain.ConvexPolygon(np.stack([[1.,0.],[0.,1.],[-1.,-1.]],axis=1))
An affine transformation can be applied to the domain, with parameters:
The direct mapping takes the form $$ x \mapsto A (x-x_0)+x_0 + v $$
def rot(t): c,s = np.cos(t),np.sin(t); return np.array(((c,-s),(s,c)))
aff = Domain.AffineTransform(cup,0.7*rot(np.pi/3.),shift=[-0.2,-0.2])
doms = (ball,box,cap,compl,cup,band,triangle,aff)
domNames = ("ball","box","cap","compl","cup","band","triangle","affine")
Let us display the domains.
aX=np.linspace(-1.2,1.2)
X = np.array(np.meshgrid(aX,aX,indexing='ij'))
h=aX[1]-aX[0]
plt.figure(figsize=(16,8))
for i,(dom,name) in enumerate(zip(doms,domNames)):
plt.subplot(2,4,1+i)
plt.contourf(*X,dom.contains(X))
plt.title(name)
plt.axis('equal')
It is possible to select the points around which the domain contains a ball of a given radius $h$, possibly negative. This predicate is only approximate.
plt.figure(figsize=(16,8))
for i,(dom,name) in enumerate(zip(doms,domNames)):
plt.subplot(2,4,1+i)
interior = dom.contains(X)
erosion = dom.contains_ball(X,2*h)
boundary_layer = np.logical_and(interior,np.logical_not(erosion))
plt.contourf(*X,boundary_layer)
plt.title(name)
plt.axis('equal')
Each domain comes equipped with a level set function, that is negative inside the domain, and positive outside. It is also guaranteed to be $1$-Lipschitz.
plt.figure(figsize=(16,8))
for i,(dom,name) in enumerate(zip(doms,domNames)):
plt.subplot(2,4,1+i)
plt.contourf(*X,dom.level(X))
plt.title(name)
plt.axis('equal')
This level set function is in general different from the signed Euclidean distance to the boundary.
If you do need the Euclidean distance function, you may consider solving an eikonal equation, or using an adequate computational geometry package.
For the design of finite difference schemes, it is important to know the distance from a given point to the domain boundary in a given direction. This is referred to as the "free way" from $x$ in the direction $v$.
#Domain = importlib.reload(Domain)
#ball = Domain.Ball(np.array([0,0]))
#box = Domain.Box([[0,1],[0,1]])
#cap = Domain.Intersection( (ball,box) )
#abox = Domain.AbsoluteComplement(box)
#aball = Domain.AbsoluteComplement(ball)
#compl = Domain.Complement(ball,box)
#cup = Domain.Union((ball,box))
#acup = Domain.Intersection((aball,abox))
#band = Domain.Band([1.,2.],[-2.,1.]) # Direction, bounds
v=np.array([1,-0.5])
plt.figure(figsize=(16,8))
for i,(dom,name) in enumerate(zip(doms,domNames)):
plt.subplot(2,4,i+1)
plt.title(name)
plt.axis('equal')
fw = dom.freeway(X,v)
if np.all(fw==np.inf): continue #Warning triggered otherwise
plt.contourf(*X,fw)
v=np.array([-1,-0.5])
plt.figure(figsize=(16,8))
for i,(dom,name) in enumerate(zip(doms,domNames)):
plt.subplot(2,4,i+1)
plt.contourf(*X,dom.freeway(X,v))
v=np.array([1,0.])
plt.figure(figsize=(16,8))
for i,(dom,name) in enumerate(zip(doms,domNames)):
plt.subplot(2,4,i+1)
plt.contourf(*X,dom.freeway(X,v))
Only Dirichlet boundary conditions are implemented at present. When a grid point falls outside the domain, the provided boundary data is used instead, as described below.
Denote by $\Omega$ the domain, and $\Omega_h$ its intersection with the cartesian grid $h Z^d$. Let $x\in \Omega_h$ be a point of the domain, and let $e \in Z^d\setminus \{0\}$ be an offset.
Define the dirichlet data, based on a domain,a function defined on (at least) the boundary, and the cartesian grid.
def bc_value(x): return x[0]+2*x[1]
bc_domain = cup
bc = Domain.Dirichlet(bc_domain,bc_value,X)
The gridscale of the domain is automatically extracted. Note that we only support gridscales which are axis independent and position independent.
bc.gridscale
0.048979591836734615
A boolean mask of interior points is constructed.
plt.title('domain interior'); plt.axis('equal')
plt.contourf(*X,bc.interior);
Actually, this is a very slightly eroded version of the domain, for numerical stability.
We choose a linear function, in order to confirm that the upwind scheme is exact in this case.
u = bc.value(X)
du = bc.DiffUpwind(u,(1,0))
du[bc.interior].max(),du[bc.interior].min()
(1.0000000000000142, 0.9999999999999833)
du,h = bc.DiffUpwind(u,(1,0),reth=True)
print("Largest and smallest h in finite differences:",np.max(h),np.min(h))
plt.title("Value of h in finite difference"); plt.axis('equal')
plt.contourf(*X,h); plt.colorbar();
Largest and smallest h in finite differences: 0.048979591836734615 0.00047605827449359595
Differentiating along the horizontal and vertical directions simultaneously.
du = bc.DiffUpwind(u,np.eye(2).astype(int))
du[0,bc.interior].max(),du[0,bc.interior].min()
(1.0000000000000142, 0.9999999999999833)
du[1,bc.interior].max(),du[1,bc.interior].min()
(2.000000000000028, 1.9999999999999858)
Again, we choose a linear function, in order to confirm that the upwind scheme is exact in this case.
du = bc.DiffCentered(u,(1,0))
du[bc.interior].max(),du[bc.interior].min()
(1.000000000000007, 0.9999999999999833)
du = bc.DiffCentered(u,np.eye(2).astype(int))
du[0,bc.interior].max(),du[0,bc.interior].min()
(1.000000000000007, 0.9999999999999833)
du[1,bc.interior].max(),du[1,bc.interior].min()
(2.000000000000028, 1.99999999999999)
The second order differences of a linear function identically vanish.
d2u = bc.Diff2(u,(1,0))
d2u[bc.interior].max(),d2u[bc.interior].min()
(1.1536608264867289e-11, -5.983034965740204e-13)
d2u = bc.Diff2(u,np.eye(2).astype(int))
d2u[0,bc.interior].max(),d2u[0,bc.interior].min()
(1.1536608264867289e-11, -5.983034965740204e-13)
d2u[1,bc.interior].max(),d2u[1,bc.interior].min()
(1.9440713794339614e-11, -5.802765675374161e-13)
We need to use quadratic polynomials to illustrate consistency in a non-trivial manner.
def bc2_value(x): return x[0]**2+x[0]*x[1]
bc2 = Domain.Dirichlet(bc_domain,bc2_value,X)
u2 = bc2.value(X)
d2u = bc2.Diff2(u2,(1,0))
d2u[bc.interior].max(),d2u[bc.interior].min()
(2.0000000000041758, 1.9999999999877116)
d2u = bc2.Diff2(u2,np.eye(2).astype(int))
d2u[0,bc.interior].max(),d2u[0,bc.interior].min()
(2.0000000000041758, 1.9999999999877116)
d2u[1,bc.interior].max(),d2u[1,bc.interior].min()
(4.265283626933634e-12, -8.961585388799509e-12)
There is an alternative way of handling boundary conditions, which is simpler but also much more crude and inaccurate. In that approach, the boundary data is provided in a neighborhood of the boundary, and the standard finite difference schemes is used.
Be warned that this alternative is both:
In view of these caveats, we refer to it as mock boundary conditions.
grid_values = bc_value(X) # Boundary conditions are evaluated in the whole domain
grid_values[ np.logical_not(bc_domain.contains(X)) ] = np.nan # Define the domain
bc_mock = Domain.MockDirichlet(grid_values,h)
In that approach, the domain is actually not taken into account when computing the finite differences.
However, the provided placeholder values are those of the boundary condition, instead of an arbitrary default. This makes sense in view of the numerical schemes implementation, see e.g. link
bc_mock.grid_values is grid_values # Returns the extended boundary conditions
True
bc.grid_values # Arbitrary placeholder
0.0