Skip to content

Fix method _test_inside_point in global smooth boundary #3

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
fhindenlang wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Fix method _test_inside_point in global smooth boundary #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 42 additions & 6 deletions pybie2d/boundaries/global_smooth_boundary/global_smooth_boundary.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
import numpy as np
import scipy as sp
import matplotlib as mpl
import scipy.signal
import warnings
import os
Expand Down Expand Up @@ -74,6 +75,7 @@ def __init__(self, x=None, y=None, c=None):
self.max_h = np.max(self.weights)
self.area = self.dt*np.sum(self.x*self.cp.imag)
self.perimeter = self.dt*np.sum(self.speed)
self.poly= mpl.path.Path(np.column_stack([self.x, self.y]))

self.defined_modules = [ 'Laplace_SLP_Self_Kress',
'Stokes_SLP_Self_Kress',
Expand Down Expand Up @@ -193,15 +195,49 @@ def Get_Close_Corrector(self, kernel, *args, **kwargs):
#########################
#### Private Methods ####
#########################

def _test_inside_point(self, candidate, eps=1e-10):

# THIS ROUTINE IS NOT RELIABLE
#def _test_inside_point(self, candidate, eps=1e-10):
# """
# Test whether the provided or generated inside point is acceptable
# returns True if the point is okay, False if its not
# """
# test_value = np.sum(self.complex_weights/(self.c-candidate))
# return np.abs(test_value - 2.0j*np.pi) < eps
## end _test_inside_point function

# THIS ROUTINE IS RELIABLE AND INDEPENDENT OF MATPLOTLIB
#def _test_inside_point(self, candidate, eps=1e-12):
# """
# Test whether the provided or generated inside point is acceptable
# returns True if the point is okay, False if its not
# Algorithm: point in polygon (the points in self.c are interpreted as closed polygon),
# the sum of angles between the vectors from the point (candidate) to two consecutive points on the boundary #
# polygon. the sum is equal 2pi only if the point lies inside. Safe way to compute the angle is using arctan2 of #(a x b)/(a.b).
# """
# a=self.c-candidate
# if (np.any(np.abs(a)<eps)):
# return False # closer than eps to a boundary point, so labelled not inside
# # angle between two "vectors" expressed as a complex number: a[i]=c[i]-cpoint,b=a[i+1]
# b=np.roll(a,-1,axis=0)
# axb=(a.real*b.imag-a.imag*b.real)
# adotb=a.real*b.real+a.imag*b.imag
# # compute sum of angles
# test_value=np.sum(np.arctan2(axb,adotb))
# return np.abs(test_value - 2.0*np.pi) < eps
## end _test_inside_point function

# THIS ROUTINE IS RELIABLE, DEPENDS OF MATPLOTLIB (same use as in `find_interior_points`)
def _test_inside_point(self, candidate):
"""
Test whether the provided or generated inside point is acceptable
returns True if the point is okay, False if its not
returns True if the point (=candidate) is okay, False if its not
Algorithm:
approximate the boundary points as a closed polygon, which is
initialized in self.poly using matplotlib.path. Then use the contains_point method
to check if the point lies inside the polygon (analogous to `find_interior_points`)
"""
test_value = np.sum(self.complex_weights/(self.c-candidate))
return np.abs(test_value - 2.0j*np.pi) < eps
# end _test_inside_point function
return self.poly.contains_point([candidate.real,candidate.imag])

def arr_check(x):
if type(x) != np.ndarray:
Expand Down
105 changes: 105 additions & 0 deletions tests/test_global_smooth_boundary.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,105 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%reload_ext autoreload\n",
"%autoreload 2\n",
"import numpy as np\n",
"import pybie2d\n",
"import matplotlib as mpl\n",
"import matplotlib.pyplot as plt\n",
"\n",
"star = pybie2d.misc.curve_descriptions.star\n",
"GSB = pybie2d.boundaries.global_smooth_boundary.global_smooth_boundary.Global_Smooth_Boundary\n",
"\n",
"def check_GSB_test_inside_points(NBfactor):\n",
" '''\n",
" Test the boundary._test_inside_point method, using points that are known to be inside/outside a star-shaped closed globally-smooth boundary\n",
" '''\n",
" \n",
" NB0=3*5\n",
" NB=NB0*NBfactor\n",
" boundary = GSB(c=((star(NB,x=1,y=1,r=1,a=0.6,f=7))))\n",
" \n",
" print(f\"NB={NB}\")\n",
" #two points on the boundary\n",
" pc = boundary.c[[1*NB//5,2*NB//3]]\n",
" eps=1e-10\n",
" t=np.linspace(eps,1-eps,30)\n",
" inner_pts=pc[0]+(pc[1]-pc[0])*t-0.5*np.sin(np.pi*t)*(pc[1].real-pc[0].real)\n",
" outer_pts=np.concatenate((pc[1]+(pc[1]-pc[0])*t+0.5*np.sin(np.pi*t)*(pc[1].real-pc[0].real),\n",
" pc[0]-(pc[1]-pc[0])*t+0.5*np.sin(np.pi*t)*(pc[1].real-pc[0].real)))\n",
"\n",
"\n",
" \n",
" check_inpts=np.array([False]*len(inner_pts))\n",
" for i,inpt in enumerate(inner_pts):\n",
" check_inpts[i]=boundary._test_inside_point(inpt)\n",
" \n",
" if not np.all(check_inpts): \n",
" print(f'test_inside_point failed for {np.sum((~check_inpts)*1)} /{30} points that are inside the boundary')\n",
" failed_inpts=inner_pts[~check_inpts]\n",
" success_inpts=inner_pts[check_inpts]\n",
"\n",
" \n",
" check_outpts=np.array([False]*len(outer_pts)) \n",
" \n",
" for i,outpt in enumerate(outer_pts):\n",
" check_outpts[i]=~boundary._test_inside_point(outpt)\n",
" \n",
" if not np.all(check_outpts): print(f'test_inside_point failed for {np.sum((~check_outpts)*1)} /{60} points that are outside the boundary')\n",
"\n",
" failed_outpts=outer_pts[~check_outpts]\n",
" success_outpts=outer_pts[check_outpts]\n",
" failed_pts=np.concatenate((failed_inpts,failed_outpts))\n",
" success_pts=np.concatenate((success_inpts,success_outpts))\n",
" return boundary.c,failed_pts,success_pts\n",
"\n",
"\n",
"\n",
"\n",
"pts,failed_pts,success_pts=check_GSB_test_inside_points(10)\n",
"fig, ax = plt.subplots(1,1,figsize=(12,12))\n",
"ax.set_aspect('equal')\n",
"\n",
"ax.plot(pts.real,pts.imag,'k-')\n",
"ax.plot(failed_pts.real,failed_pts.imag,'rx')\n",
"ax.plot(success_pts.real,success_pts.imag,'bx') \n",
"\n",
"assert(len(failed_pts)==0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"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": 2
}