agd.Domain

  1# Copyright 2020 Jean-Marie Mirebeau, University Paris-Sud, CNRS, University Paris-Saclay
  2# Distributed WITHOUT ANY WARRANTY. Licensed under the Apache License, Version 2.0, see http://www.apache.org/licenses/LICENSE-2.0
  3
  4"""
  5This module allows to define domains of $R^d$ by combinations of elementary shapes, 
  6and to compute finite differences within these domains with Dirichlet boundary conditions.
  7"""
  8
  9from . import FiniteDifferences as fd
 10from . import AutomaticDifferentiation as ad
 11from . import LinearParallel as lp
 12
 13import numpy as np
 14from functools import reduce
 15from itertools import cycle
 16
 17class Domain:
 18	"""
 19	This class represents a domain from which one can query 
 20	a level set function, the boundary distance in a given direction,
 21	and some related methods.
 22	"""
 23
 24	def __init__(self):
 25		pass
 26
 27	def level(self,x):
 28		"""
 29		A level set function, negative inside the domain, positive outside.
 30		Guaranteed to be 1-Lipschitz.
 31		"""
 32		raise ValueError("""Domain level set function must be specialized""")
 33
 34	def contains(self,x):
 35		"""
 36		Wether x lies inside the domain.
 37		"""
 38		return self.level(x)<0
 39
 40	def intervals(self,x,v):
 41		r"""
 42		A union of disjoint intervals, sorted in increasing order, 
 43		$$
 44			] a_0,b_0 [ \cup ] a_1,b_1 [ \cup ... \cup ] a_{n-1},b_{n-1} [
 45		$$
 46		such that $x+t v$ lies in the domain iff $t$ lies on one of these intervals.
 47		"""
 48		raise ValueError("""Domain intervals function must be specialized""")
 49
 50	def freeway(self,x,v):
 51		r"""
 52		Output : Least $t\geq 0$ such that $x+tv$ intersects the boundary.
 53		"""
 54		a,b = self.intervals(x,v)
 55		a[a<0]=np.inf
 56		b[b<0]=np.inf
 57		return np.minimum(a.min(axis=0),b.min(axis=0))
 58
 59	def contains_ball(self,x,h,ball_pattern=None):
 60		"""
 61		Wether the domain contains the ball of center $x$ and radius $h$.
 62		Approximate predicate based on sampling, using ball_pattern.
 63		"""
 64		if h==0.: 	return self.contains(x)
 65		if h<0.:	return AbsoluteComplement(self).contains_ball(x,-h)
 66		level = self.level(x)
 67		inside = level<0
 68
 69		# Fix boundary layer
 70		mask = np.logical_and(inside,level>-h) # Recall level is 1-Lipschitz
 71		xm=x[:,mask]
 72		if ball_pattern is None: ball_pattern = self.ball_pattern()
 73		ball_pattern = fd.as_field(ball_pattern,xm.shape[1:],conditional=False)
 74		xb = np.expand_dims(xm,axis=1)+h*ball_pattern
 75		inside[mask] = np.all(self.contains(xb),axis=0)
 76		return inside
 77
 78	def ball_pattern(self,vdim=None,r=1):
 79		"""
 80		Produces a sampling pattern in a ball, to be used in the contains_ball predicate.
 81		"""
 82		if vdim is None: vdim = self.vdim
 83		if vdim is None: raise ValueError("Unspecified dimension")
 84		aX = list(range(-r,r+1))
 85		X = np.asarray(np.meshgrid(*(aX,)*vdim,indexing='ij')).reshape(vdim,-1)
 86		nX = np.sqrt((X**2).sum(axis=0))
 87		pos = nX>0
 88		return X[:,pos]/nX[pos]	
 89
 90	@property
 91	def vdim(self):
 92		"""Dimension of the embedding space"""
 93		return None
 94	
 95
 96class WholeSpace(Domain):
 97	"""
 98	This class represents the full space $R^d$.
 99	"""
100	def level(self,x):
101		return np.full(x.shape[1:],-np.inf)
102
103	def intervals(self,x,v):
104		shape = (1,)+x.shape[1:]
105		return np.full(shape,-np.inf),np.full(shape,np.inf)
106
107class Ball(Domain):
108	"""
109	This class represents a ball shaped domain.
110
111	__init__ arguments : 
112	- center (optional), array : the center of the ball.
113	- radius (optional), scalar : the radius of the ball.
114	Defaults to defining the unit two-dimensional ball.
115	"""
116
117	def __init__(self,center=(0.,0.),radius=1.):
118		super(Ball,self).__init__()
119		self.center = ad.asarray(center)
120		self.radius = radius
121
122	@property
123	def vdim(self): return len(self.center)
124
125	def _centered(self,x):
126		_center = fd.as_field(self.center,x.shape[1:],conditional=False)
127		return x-_center
128
129	def level(self,x):
130		_x = self._centered(x)
131		return ad.Optimization.norm(_x,ord=2,axis=0)-self.radius
132
133	def intervals(self,x,v):
134		if v.shape!=x.shape: v=fd.as_field(v,x.shape[1:],conditional=False)
135		xc = self._centered(x)
136
137		begin = np.full(x.shape[1:],np.inf)
138		end = begin.copy()
139
140		# Solve |x+hv|^2=r^2, which is a quadratic equation a h^2 + 2 b h + c =0
141		a = lp.dot_VV(v,v)
142		b = lp.dot_VV(xc,v)
143		c = lp.dot_VV(xc,xc)-self.radius**2
144
145		delta = b*b-a*c
146
147		pos = np.logical_and(a>0,delta>0)
148		a,b,delta = (e[pos] for e in (a,b,delta))
149
150		sdelta = np.sqrt(delta)
151		begin[pos] = (-b-sdelta)/a
152		end[pos] = (-b+sdelta)/a
153
154		begin,end = (np.expand_dims(e,axis=0) for e in (begin,end))
155		return begin,end
156
157class Box(Domain):
158	r"""
159	This class represents a box shaped domain in $R^k$, mathematically defined 
160	as the product
161	$$
162		[a_1,b_1] \times ... \times [a_k,b_k]
163	$$
164	of some intervals. 
165
166
167	__init__ argument : 
168	- sides (optional) : [[a1,b1], ..., [ak,bk]] the intervals defining the box.
169	 Defaults to defining the two dimensional unit square.
170	"""
171
172	def __init__(self,sides = ((0.,1.),(0.,1.)) ):
173		super(Box,self).__init__()
174		sides = ad.asarray(sides)
175		self._sides = sides
176		self._center = sides.sum(axis=1)/2.
177		self._hlen = (sides[:,1]-sides[:,0])/2.
178	
179	@property
180	def sides(self): 
181		"""
182		The intervals $[[a_1,b_1], ..., [a_k,b_k]]$ defining the box.
183		"""
184		return self._sides
185
186	@property
187	def center(self): 
188		"""
189		The center $[ (a_1+b_1)/2, ..., (a_k+b_k)/2 ]$ of the box.
190		"""
191		return self._center
192
193	@property
194	def edgelengths(self): 
195		"""
196		The edge lengths $[ b_1-a_1, ..., b_k-a_k ]$ of the box.
197		"""
198		return 2.*self._hlen
199	@property
200	def vdim(self): return len(self.center)
201
202	def _centered(self,x,signs=False):
203		center = fd.as_field(self.center,x.shape[1:],conditional=False)
204		xc = x-center
205		axc = np.abs(xc)
206		if not signs: return axc
207		sxc = 2*(xc>=0)-1 # Non-vanishing sign 
208		return axc,sxc
209
210	def level(self,x):
211		hlen = fd.as_field(self._hlen,x.shape[1:],conditional=False)
212		return (self._centered(x) - hlen).max(axis=0)
213
214	def intervals(self,x,v):
215		shape = x.shape[1:]
216		if v.shape!=x.shape: v=fd.as_field(v,shape,conditional=False)
217		hlen = fd.as_field(self._hlen,shape,conditional=False)
218		xc,signs = self._centered(x,signs=True)
219		vc = v*signs
220
221		a=np.full(x.shape,-np.inf)
222		b=np.full(x.shape, np.inf)
223
224		# Compute the interval corresponding to each axis
225		pos = vc!=0
226		hlen_,xc_,vc_ = (e[pos] for e in (hlen,xc,vc))
227
228		a[pos] = (-hlen_-xc_)/vc_
229		b[pos] = ( hlen_-xc_)/vc_
230
231		# Deal with offsets parallel to axes
232		pos = np.logical_and(vc==0.,np.abs(xc)>hlen)
233		a[pos] = np.inf 
234
235		# Intersect intervals corresponding to different axes
236
237		a,b = np.minimum(a,b),np.maximum(a,b)
238		a,b = a.max(axis=0),b.min(axis=0)
239		a,b = (ad.asarray(e) for e in (a,b)) 
240
241		# Normalize empty intervals
242		pos = a>b
243		a[pos]=np.inf
244		b[pos]=np.inf
245		a,b = (np.expand_dims(e,axis=0) for e in (a,b))
246		return a,b
247
248class AbsoluteComplement(Domain):
249	r"""
250	This class represents the complement $R^d \setminus \Omega$, in the entire space,
251	 of a given domain $\Omega\subset R^d$.
252
253	__init__ argument: 
254	- dom : Domain of which to take the complement.
255	"""
256	def __init__(self,dom):
257		super(AbsoluteComplement,self).__init__()
258		self.dom = dom
259
260	@property
261	def vdim(self): return self.dom.vdim
262	def contains(self,x):
263		return np.logical_not(self.dom.contains(x))
264	def level(self,x):
265		return -self.dom.level(x)
266	def freeway(self,x,v):
267		return self.dom.freeway(x,v)
268	def intervals(self,x,v):
269		a,b = self.dom.intervals(x,v)
270		inf = np.full((1,)+x.shape[1:],np.inf)
271		return np.concatenate((-inf,b),axis=0),np.concatenate((a,inf),axis=0)
272
273class Intersection(Domain):
274	"""
275	This class represents an intersection of several subdomains.
276
277	__init__ arguments : 
278	- *doms : domains to intersect.
279	"""
280	def __init__(self,*doms):
281		super(Intersection,self).__init__()
282		self.doms=doms
283
284	@property
285	def vdim(self): 
286		vdims=[dom.vdim for dom in self.doms if dom.vdim is not None]
287		if len(vdims)==0: return None
288		vdim = vdims[0]
289		assert vdims.count(vdim)==len(vdims)
290		return vdim
291	
292	def contains(self,x):
293		containss = [dom.contains(x) for dom in self.doms]
294		return reduce(np.logical_and,containss)
295
296	def level(self,x):
297		levels = [dom.level(x) for dom in self.doms]
298		return reduce(np.maximum,levels)
299
300	def intervals(self,x,v):
301		intervalss = [dom.intervals(x,v) for dom in self.doms]
302		begs = [a for a,b in intervalss]
303		ends = [b for a,b in intervalss]
304
305		# Shortcut for the convex case, quite common
306		if all(len(a)==1 for a,b in intervalss):
307			beg = reduce(np.maximum,begs)
308			end = reduce(np.minimum,ends)
309			pos = beg>end
310			beg[pos]=np.inf
311			end[pos]=np.inf
312			return beg,end
313
314		# General non-convex case
315		shape = x.shape[1:]
316
317		# Prepare the result
318		begr = []
319		endr = []
320		val = np.full(shape,-np.inf)
321		inds = np.full( (len(self.doms),)+x.shape[1:], 0)
322		
323		def tax(arr,ind,valid):
324			return np.squeeze(np.take_along_axis(arr[:,valid],np.expand_dims(ind[valid],axis=0),axis=0),axis=0)
325
326		counter=0
327		while True:
328			# Find the next beginning
329			unchanged=0
330			for it,(beg,end) in cycle(enumerate(zip(begs,ends))):
331				ind=ad.asarray(inds[it])
332				valid = ind<len(end)
333				pos = np.full(shape,False)
334				endiv = tax(end,ind,valid)
335				pos[valid] = np.logical_and(val[valid]>=endiv,endiv!=np.inf)
336
337				if pos.any():	
338					unchanged=0
339				else:			
340					unchanged+=1
341					if unchanged>=len(begs):
342						break
343
344				inds[it,pos]+=1; ind=inds[it] 
345				valid = ind<len(end)
346				val[valid] = np.maximum(val[valid],tax(beg,ind,valid))
347				val[np.logical_not(valid)]=np.inf
348
349				counter+=1
350				assert(counter<=100)
351
352			#Exit if nobody's valid
353			if (val==np.inf).all():
354				break
355			begr.append(val)
356
357			# Find the next end
358			val=np.full(shape,np.inf)
359			for end,ind in zip(ends,inds):
360				valid = ind<len(end)
361				val[valid] = np.minimum(val[valid],tax(end,ind,valid))
362
363			endr.append(val.copy())
364
365		return np.asarray(begr),np.asarray(endr)
366	
367def Complement(dom1,dom2):
368	r"""
369	Relative complement $\Omega_1 \setminus \Omega_2$ of two given domains 
370	$\Omega_1,\Omega_2 \subset R^d$.
371	"""
372	return Intersection(dom1,AbsoluteComplement(dom2))
373
374def Union(*doms):
375	"""
376	Union of several domains.
377	"""
378	return AbsoluteComplement(Intersection(*[AbsoluteComplement(dom) for dom in doms]))
379
380class Band(Domain):
381	r"""
382	Defines a banded domain in space, in between two parallel hyperplanes.
383	$$
384		l_b < < x,v> < u_b,
385	$$
386	where $v \in R^d$ is a direction, and $l_b,u_b$ are a lower bound and an upper bound.
387
388	__init__ arguments : 
389	- direction : array of shape $(d,)$, where $d$ is the ambient space dimension.
390	- bounds : $[l_b, u_b]$, the lower bound and upper bound.
391	"""
392
393	def __init__(self,direction,bounds):
394		super(Band,self).__init__()
395		self.direction,self.bounds = (ad.asarray(e) for e in (direction,bounds))
396		norm = ad.Optimization.norm(self.direction,ord=2)
397		if norm!=0.:
398			self.direction/=norm
399			self.bounds/=norm
400
401	@property
402	def vdim(self): return len(self.direction)
403	
404	def _dotdir(self,x):
405		direction = fd.as_field(self.direction,x.shape[1:],conditional=False)
406		return lp.dot_VV(x,direction)
407
408	def level(self,x):
409		xd = self._dotdir(x)
410		return np.maximum(self.bounds[0]-xd,xd-self.bounds[1])
411
412	def intervals(self,x,v):
413		xd = self._dotdir(x)
414		vd = self._dotdir(v)
415		vd = fd.as_field(vd,xd.shape)
416		a = np.full(xd.shape,-np.inf)
417		b = np.full(xd.shape, np.inf)
418
419		# Non degenerate case
420		mask = vd!=0
421		xdm,vdm = xd[mask],vd[mask]
422		a[mask] = (self.bounds[0]-xdm)/vdm
423		b[mask] = (self.bounds[1]-xdm)/vdm
424		# Handle the case where vd=0
425		inside = np.logical_and(self.bounds[0]<xd,xd<self.bounds[1])
426		a[np.logical_and(vd==0,np.logical_not(inside))]=np.inf
427
428		a,b = np.minimum(a,b),np.maximum(a,b)
429		a,b = (np.expand_dims(e,axis=0) for e in (a,b))
430		return a,b
431
432def ConvexPolygon(pts):
433	"""
434	Defines a two dimensional *convex* polygonal domain from its vertices.
435
436	__init__ arguments :
437	- pts : array of shape (2,n). The polygon vertices, given in trigonometric order.
438	"""
439	def params(p,q):
440		pq = q-p
441		direction = np.asarray([pq[1],-pq[0]])
442		lower_bound = np.dot(direction,p)
443		return direction,[lower_bound,np.inf]
444
445	pts = ad.asarray(pts)
446	assert len(pts)==2
447	return Intersection(*[Band(*params(p,q)) for p,q in zip(pts.T,np.roll(pts,1,axis=1).T)])
448
449class AffineTransform(Domain):
450	"""
451	Defines a domain which is the image of another domain
452	by the affine transformation, 
453	$$
454		x' = A x + b
455	$$
456
457	Inputs : 
458	- `mult` (optional) : the multiplier $A$, which is either a scalar or a matrix 
459	of shape $(d,d)$. (Defaults to the identity matrix.)
460	- `shift` (optional) : the vector $b$. (Defaults to zero.)
461	- `center` (optional) : reference point for the transformation. (Defaults to zero.)
462	"""
463
464	def __init__(self,dom,mult=None,shift=None,center=None):
465		super(AffineTransform,self).__init__()
466		self.dom=dom
467		if mult is not None: mult = ad.asarray(mult)
468		self._mult = mult
469
470		if shift is not None: shift = ad.asarray(shift)
471		if center is not None:
472			center=ad.asarray(center)
473			shift2=center-self.forward(center,linear=True)
474			shift = shift2 if shift is None else shift+shift2
475
476		self._shift = shift
477		self._mult_inv = (None if mult is None else 
478			(ad.asarray(1./mult) if mult.ndim==0 else np.linalg.inv(mult) ) )
479		self._mult_inv_norm = (1. if self._mult_inv is None 
480			else np.linalg.norm(self._mult_inv,ord=2) )
481
482	@property
483	def vdim(self): return self.dom.vdim
484
485	def forward(self,x,linear=False):
486		"""
487		Forward affine transformation, from the original domain to the transformed one.
488		"""
489		x=x.copy()
490		shape = x.shape[1:]
491
492		mult = self._mult
493		if mult is None: 	pass
494		elif mult.ndim==0:	x*=mult
495		else: x=lp.dot_AV(fd.as_field(mult,shape,conditional=False),x)
496
497		if not linear and self._shift is not None:
498			x+=fd.as_field(shift,shape,conditional=False)
499		
500		return x
501
502	def reverse(self,x,linear=False):
503		"""
504		Reverse affine transformation, from the transformed domain to the original one.
505		"""
506		x=x.copy()
507		shape = x.shape[1:]
508
509		shift = self._shift
510		if (shift is None) or linear: 	pass
511		else:	x-=fd.as_field(shift,shape,conditional=False)
512
513		mult = self._mult_inv
514		if mult is None: 	pass
515		elif mult.ndim==0:	x*=mult
516		else: x = lp.dot_AV(fd.as_field(mult,shape,conditional=False),x)
517
518		return x
519
520	def contains(self,x):
521		return self.dom.contains(self.reverse(x))
522	def level(self,x):
523		return self.dom.level(self.reverse(x))/self._mult_inv_norm
524	def intervals(self,x,v):
525		return self.dom.intervals(self.reverse(x),self.reverse(v,linear=True))
526	def freeway(self,x,v):
527		return self.dom.freeway(self.reverse(x),self.reverse(v,linear=True))
528
529class Dirichlet:
530	"""
531	Implements Dirichlet boundary conditions.
532	When computing finite differences, values queried outside the domain interior
533	are replaced with values obtained on the boundary along the given direction.
534
535	__init__ arguments :
536	- domain: geometrical description of the domain 
537	- value: a scalar or map yielding the value of the boundary conditions
538	- grid: the cartesian grid. Ex: np.array(np.meshgrid(aX,aY,indexing='ij')) for suitable aX,aY
539
540	- interior (optional): the points regarded as interior to the domain. 
541	- interior_radius (optional): sets
542		interior = domain.contains_ball(interior_radius)
543
544	- grid_values (optional): placeholder values to be used on the grid.
545		Either an array of values, or a function
546
547	member fields : 
548	- interior : mask for the discretization grid points inside the domain.
549	"""
550
551	def __init__(self,domain,value,grid,
552		interior_radius=None,interior=None,grid_values=0.):
553		self.domain = domain
554
555		if isinstance(value,float):
556			self.value = lambda x:value
557		else:
558			self.value = value
559
560		self.grid = ad.asarray(grid)
561		self.gridscale = self._gridscale(self.grid)
562
563		if interior is not None:
564			self.interior = interior
565		else:
566			if interior_radius is None:
567#				interior_radius = 0.5 * self.gridscale
568				# Tiny value, for numerical stability. Ideally 0.
569				interior_radius = self.gridscale * 1e-8 
570			self.interior = self.domain.contains_ball(self.grid,interior_radius)
571
572		if isinstance(grid_values,float) or ad.isndarray(grid_values):
573			self.grid_values = grid_values
574		else:
575			self.grid_values = grid_values(self.grid)
576
577	
578	@staticmethod
579	def _gridscale(grid):
580		dim = len(grid)
581		assert(grid.ndim==1+dim)
582		x0 = grid.__getitem__((slice(None),)+(0,)*dim)
583		x1 = grid.__getitem__((slice(None),)+(1,)*dim)
584		delta = np.abs(x1-x0)
585		hmin,hmax = np.min(delta),np.max(delta)
586		if hmax>=hmin*(1+1e-8):
587			raise ValueError("Error : gridscale is axis dependent")
588		return hmax
589
590	@property 
591	def shape(self): 
592		"""
593		The shape of the domain (discretization grid).
594		"""
595		return self.grid.shape[1:]
596
597	@property
598	def vdim(self): 
599		"""
600		The dimension of the vector space containing the domain.
601		"""
602		return len(self.grid)	
603	def as_field(self,arr,**kwargs): 
604		"""
605		Adds trailing dimensions, and broadcasts arr, if necessary,
606		so that the shape ends with the domain shape. 
607		"""
608		return fd.as_field(arr,self.shape,**kwargs)
609
610	@property
611	def not_interior(self): 
612		"""
613		Mask for the grid points outside the domain.
614		"""
615		return np.logical_not(self.interior)
616	
617	@property
618	def Mock(self):
619		"""
620		Returns mock Dirichlet boundary conditions obtained by evaluating 
621		the boundary condition outside the interior.
622		"""
623		bc_values = np.full(self.grid.shape[1:],np.nan)
624		bc_values[self.not_interior] = self.value(self.grid[:,self.not_interior])
625		return MockDirichlet(bc_values,self.gridscale)
626
627	@property
628	def _ExteriorNaNs(self):
629		result = np.zeros(self.grid.shape[1:])
630		result[self.not_interior]=np.nan
631		return result
632
633	def _BoundaryLayer(self,u,du):
634		"""
635		Returns positions at which u is defined but du is not.
636		"""
637		return np.logical_and.reduce([
638			np.isnan(du),
639			np.broadcast_to(self.interior,du.shape),
640			np.broadcast_to(np.logical_not(np.isnan(u)),du.shape)
641			])
642
643	def _DiffUpwindDirichlet(self,u,offsets,grid,reth):
644		"""
645		Returns first order finite differences w.r.t. boundary value.
646		"""
647		h = self.domain.freeway(grid,offsets)
648		x = grid+h*offsets
649
650		if np.any(np.isnan(x)):
651			bad = np.isnan(x.sum(axis=0))
652			print("bad",grid[:,bad],offsets[:,bad],h[bad],x[:,bad])
653
654		xvalues = self.value(x)
655		result = (xvalues-u)/h
656		return (result,h) if reth else result
657
658
659
660	def DiffUpwind(self,u,offsets,reth=False):
661		"""
662		First order upwind finite differences of u, along the offsets direction. 
663		Uses boundary values when needed.
664		Input : 
665		- u : array with the domain shape.
666		- offsets :  array of integers, with shape 
667		(vdim, n1,...,nk) or (vdim, n1,...,nk, shape) 
668		where vdim,shape is the ambient space dimension and domain shape, 
669		and n1,...,nk are arbitrary.
670		- reth (optional) : wether to return the grid scale used (differs from the 
671		grid scale on the neighborhood of the boundary).
672		"""
673		assert(isinstance(reth,bool))
674		grid=self.grid
675		du = fd.DiffUpwind(u+self._ExteriorNaNs,offsets,self.gridscale)
676		mask = self._BoundaryLayer(u,du)
677		offsets = fd.as_field(np.asarray(offsets),u.shape)
678		um = np.broadcast_to(u,offsets.shape[1:])[mask]
679		om = offsets[:,mask]
680		gm = np.broadcast_to(grid.reshape( 
681			(len(grid),)+(1,)*(offsets.ndim-grid.ndim)+u.shape),offsets.shape)[:,mask]
682#		um,om,gm = u[mask], offsets[:,mask], grid[:,mask]
683		if not reth: 
684			du[mask] = self._DiffUpwindDirichlet(um,om,gm,reth=reth)
685			return du
686		else: 
687			hr = np.full(offsets.shape[1:],self.gridscale)
688			du[mask],hr[mask] = self._DiffUpwindDirichlet(um,om,gm,reth=reth)
689			return du,hr
690
691	def DiffCentered(self,u,offsets):
692		"""
693		Centered finite differences of u, along the offsets direction, 
694		computed using the second order accurate centered scheme.
695		Falls back to upwind finite differences close to the boundary.
696		"""
697		du = fd.DiffCentered(u+self._ExteriorNaNs,offsets,self.gridscale)
698		mask = self._BoundaryLayer(u,du)
699		du1 = self.DiffUpwind(u,offsets)
700		du[mask]=du1[mask]
701		#du[mask] = self.DiffUpwind(u[mask],offsets[:,mask],h,grid[:,mask])
702		return du
703
704
705
706	def Diff2(self,u,offsets):
707		"""
708		Second order finite differences of u, along the offsets direction.
709		Second order accurate in the interior, 
710		but only first order accurate at the boundary.
711		"""
712		du0,h0 = self.DiffUpwind(u, offsets,		  reth=True)
713		du1,h1 = self.DiffUpwind(u,-np.asarray(offsets),reth=True)
714
715		return (du0+du1)*(2./(h0+h1))
716
717class MockDirichlet:
718	"""
719	Implements a crude version of Dirichlet boundary conditions, 
720	where the boundary conditions are given on the full domain complement.
721
722	(No geometrical computations involved.)
723
724	__init__ arguments : 
725	- grid_values : the Dirichlet boundary conditions. 
726	 Please set to np.nan the values interior to domain.
727	- gridscale : the discretization grid scale.
728	- padding : the padding values to be used outside the discretization grid.
729	"""
730
731	def __init__(self,grid_values,gridscale,padding=np.nan):
732		if isinstance(grid_values,tuple):
733			grid_values = np.full(grid_values,np.nan)
734		self.grid_values=grid_values
735		self.gridscale=gridscale
736		self.padding = padding
737
738	@property
739	def interior(self): 
740		"""
741		The discretization grid points inside the domain.
742		"""
743		return np.isnan(self.grid_values)
744	@property
745	def not_interior(self): 
746		"""
747		The discretization grid points outside the domain.
748		"""
749		return np.logical_not(self.interior)
750
751	@property
752	def vdim(self): 
753		"""
754		The dimension of the vector space containing the domain.
755		"""
756		return self.grid_values.ndim	
757
758	@property
759	def shape(self): 
760		"""
761		The shape of the domain discretization grid.
762		"""
763		return self.grid_values.shape
764
765	def as_field(self,arr,**kwargs): 
766		"""
767		Adds trailing dimensions, and broadcasts arr, if necessary,
768		so that the shape ends with the discretization grid shape. 
769		"""
770		return fd.as_field(arr,self.shape,**kwargs)
771
772	def DiffUpwind(self,u,offsets,**kwargs): 
773		"""
774		First order upwind finite differences of u, along the offsets direction. 
775		Uses grid_values outside the domain interior.
776		"""
777		return fd.DiffUpwind(u,offsets,self.gridscale,padding=self.padding,**kwargs)
778
779	def DiffCentered(self,u,offsets): 
780		"""
781		First order centered finite differences of u, along the offsets direction. 
782		Uses grid_values outside the domain interior.
783		"""
784		return fd.DiffCentered(u,offsets,self.gridscale,padding=self.padding)
785
786	def Diff2(self,u,offsets,*args): 
787		"""
788		Second order order finite differences of u, along the offsets direction. 
789		Uses grid_values outside the domain interior.
790		"""
791		return fd.Diff2(u,offsets,self.gridscale,*args,padding=self.padding)
792
793	
794	
795