diff --git a/src/yapsut/__init__.py b/src/yapsut/__init__.py index 32a5a29594a8a48ed620ddc11063b193e9a9b0fd..58fe84160441ac74c171cea05aad0ca1a8f8cd6a 100644 --- a/src/yapsut/__init__.py +++ b/src/yapsut/__init__.py @@ -11,6 +11,7 @@ from .profiling import timer,timerDict from .struct import struct from .DummyCLI import DummyCLI from .nearly_even_split import nearly_even_split +from .geometry import ElipsoidTriangulation #from import .grep #import .intervalls diff --git a/src/yapsut/geometry.py b/src/yapsut/geometry.py new file mode 100644 index 0000000000000000000000000000000000000000..e439d593ce0ee1b11caa3bc853d780e210a735ce --- /dev/null +++ b/src/yapsut/geometry.py @@ -0,0 +1,280 @@ +import numpy as np +class ElipsoidTriangulation : + """ creates an elipsoid triangulated """ + @property + def axis(self) : + """ elipsoid axis""" + return self._axis + @property + def shape(self) : + """ grid shape""" + return self._shape + @property + def nphi(self) : + """ nphi shape""" + return self._nphi + @property + def ntheta(self) : + """ ntheta shape""" + return self._ntheta + @property + def phi(self) : + """ phi list""" + return self._phi + @property + def theta(self) : + """ theta list""" + return self._theta + @property + def simplices(self) : + """ list of triangles (simplices)""" + return self._Triangles + @property + def simplices(self) : + """ list of triangles (simplices)""" + return self._Triangles + @property + def polarS(self) : + """ list of polar coordinates (for each point in each simplice)""" + return self._TrAng + @property + def polarS(self) : + """ list of polar coordinates (simplices)""" + return self._TrAng + @property + def vecS(self) : + """ list of cartesian coordinates for each triangle""" + return self._TrVec + @property + def verS(self) : + """ list of cartesian coordinates for each triangle as versors""" + return self._TrVec + @property + def normalS(self) : + """ list of oriented normals for each triangle""" + return self._TrNorm + @property + def polar_normalS(self) : + """ list of intercepts for each triangle, polar coordinates""" + return self._TrNormAng + @property + def interceptS(self) : + """ list of intercepts for each triangle""" + return self._TrItcp + @property + def areaS(self) : + """ list of areas for each triangle""" + return self._TrArea + @property + def areaS(self) : + """ list of areas for each triangle""" + return self._TrArea + @property + def vecCpsS(self) : + """ list of coplanar coordinates for each triangle""" + return self._TrVecCpl + # + def __len__(self) : + return self._N + # + def __init__(self,a,b,c,nphi,ntheta) : + self._axis=np.array([a,b,c]) + self._nphi=nphi + self._ntheta=ntheta + self._shape=(ntheta,nphi) + + self._phi=np.linspace(0,1,self._nphi)*2*np.pi + self._theta=np.linspace(0,1,self._ntheta)*np.pi + + phi=self._phi + theta=self._theta + + # + # map of triangles + Triang=[] + + # north pole triangles + itheta=0 + iP=[0,0] + for iiph in np.arange(nphi-1) : + iph0=np.mod(iiph,nphi) + iph1=np.mod(iiph+1,nphi) + # + T=[iP,[iph0,itheta+1],[iph1,itheta+1],'0,1'] + Triang.append(T) + + # bands triangles + Theta=[] + for itheta in np.arange(1,ntheta-2) : + itheta1=itheta+1 + ss=str(itheta)+'-'+str(itheta+1) + for iiph in np.arange(nphi-1) : + iph0=np.mod(iiph,nphi) + iph1=np.mod(iiph+1,nphi) + # + T=[[iph0,itheta],[iph0,itheta1],[iph1,itheta1],ss] + Triang.append(T) + # + T=[[iph1,itheta1],[iph1,itheta],[iph0,itheta],ss] + Triang.append(T) + + # south pole triangles + itheta=ntheta-2 + itheta1=ntheta-1 + iP=[0,itheta1] + ss='0-'+str(itheta1) + for iiph in np.arange(nphi-1) : + iph0=np.mod(iiph,nphi) + iph1=np.mod(iiph+1,nphi) + # + T=[[iph0,itheta],iP,[iph1,itheta],ss] + Triang.append(T) + self._Triangles=Triang + self._N=len(Triang) + + # convert map of triangles in angles + TrAng=[] + for BL in Triang : + T=[] + for iT in BL[:3] : + T.append([phi[iT[0]],theta[iT[1]]]) + TrAng.append(T) + TrAng=np.array(TrAng) + self._Angles=TrAng + + # convert map of triangles in vectors + TrVec=[] + for BL in Triang : + T=[] + for iT in BL[:3] : + ph=phi[iT[0]] + cph=np.cos(ph) ; sph=np.sin(ph) + th=theta[iT[1]] + cth=np.cos(th) ; sth=np.sin(th) + T.append([cph*sth,sph*sth,c*cth]) + TrVec.append(T) + TrVec=np.array(TrVec) + self._TrVec=TrVec + + # convert map of triangles in versors + TrVers=[] + for BL in TrVec : + T=[] + for iT in BL[:3] : + T.append(iT/(iT.dot(iT))**0.5) + TrVers.append(T) + TrVers=np.array(TrVers) + self._TrVers=TrVers + + # normals of triangles as vectors, angles and intercept of planes + # TrDprod is for each triangle the n.dot(r_i) with r_i are the triangles vertex + TrDprod=[] + TrNorm=[] + TrNormAng=[] + TrItcp=[] + for BL in TrVec : + r1,r2,r3=BL + d21=r2-r1 + d31=r3-r1 + cdprod=np.cross(d21,d31) + n=cdprod/(cdprod.dot(cdprod))**0.5 + # + sgn=1 + a1=n.dot(r1)/r1.dot(r1)**0.5 + a2=n.dot(r2)/r2.dot(r2)**0.5 + a3=n.dot(r3)/r3.dot(r3)**0.5 + if a1<0 or a2<0 or a3<0 : + n=-1*n + sgn=-1 + a1=n.dot(r1)/r1.dot(r1)**0.5 + a2=n.dot(r2)/r2.dot(r2)**0.5 + a3=n.dot(r3)/r3.dot(r3)**0.5 + TrItcp.append([-n.dot(r1),sgn]) + TrNorm.append(n) + TrDprod.append([a1,a2,a3]) + # + TrNormAng.append([np.arctan2(n[1],n[0]),np.arccos(n[2])]) + TrItcp=np.array(TrItcp) + TrNorm=np.array(TrNorm) + TrDprod=np.array(TrDprod) + TrNormAng=np.array(TrNormAng) + self._TrItcp=TrItcp + self._TrNorm=TrNorm + self._TrDprod=TrDprod + self._TrNormAng=TrNormAng + + # computes the are of triangles and VecCpl, triangles vertex in coplanar representation + TrVecCpl=[] + TrArea = [] + for ibl in np.arange(len(TrNorm)) : + n=TrNorm[ibl] + trv=TrVec[ibl] + # + nlo,nla=TrNormAng[ibl] + cnlo=np.cos(nlo) + cnla=np.cos(nla) + snlo=np.sin(nlo) + snla=np.sin(nla) + # + Mlo=np.array([[cnlo,-snlo,0],[snlo,cnlo,0],[0,0,1]]) + Mla=np.array([[cnla,0,snla],[0,1,0],[-snla,0,cnla]]) + #Mlo.dot(Mla.dot(np.array([0,0,1]))) + Mlo.dot(Mla.dot(np.array([0,0,1]))) + nrot=Mla.T.dot(Mlo.T.dot(n)) + # + # method of determinant and closed form + # rotate triangle in coplanar coordinates (z = constant) + ctr=np.array([Mla.T.dot(Mlo.T.dot(k)) for k in trv]) + cr1,cr2,cr3=ctr + TrVecCpl.append(ctr) + # computes area closed form + area1=0.5*abs(cr1[0]*(cr2[1]-cr3[1])+cr2[0]*(cr3[1]-cr1[1])+cr3[0]*(cr1[1]-cr2[1])) + # computes area method of determinant + ## z forced to be 1 + #ctr[:,2]=1 + #area2=0.5*scipy.linalg.det(ctr) + #TrArea.append([area2,area1]) + TrArea.append(area1) + TrArea=np.array(TrArea).T + TrVecCpl=np.array(TrVecCpl) + self._TrArea=TrArea + self._TrVecCpl=TrVecCpl + # + def dotNormals(self,v) : + """ computes dot product between normals of simplices and a vector v """ + return np.array([v.dot(k) for k in self.normalS]) + # + # da verificare + #def rotateSimplices(self,lo,la) : + #""" computes rotations of simplices + #lo = longitude [rad] + #la = latitude [rad] + #""" + #cnlo=np.cos(lo) + #cnla=np.cos(la) + #snlo=np.sin(lo) + #snla=np.sin(la) + ## + #Mlo=np.array([[cnlo,-snlo,0],[snlo,cnlo,0],[0,0,1]]) + #Mla=np.array([[cnla,0,snla],[0,1,0],[-snla,0,cnla]]) + #out=[] + #for T in self.vecS : + #out.append(Mla.T.dot(Mlo.T.dot(k)) for k in T) + #return np.array(out) + # + # da verificare + #def rotateNormals(self,lo,la) : + #""" computes rotations of normals + #lo = longitude [rad] + #la = latitude [rad] + #""" + #cnlo=np.cos(lo) + #cnla=np.cos(la) + #snlo=np.sin(lo) + #snla=np.sin(la) + ## + #Mlo=np.array([[cnlo,-snlo,0],[snlo,cnlo,0],[0,0,1]]) + #Mla=np.array([[cnla,0,snla],[0,1,0],[-snla,0,cnla]]) + ## + #return np.array([Mla.T.dot(Mlo.T.dot(k)) for k in self.normalS]) +