compute_distance_forward<-function(p,n,mu){# for the direction mu we compute the distance always forward h<-(n[mu]-p[mu])if (h<=0){ h<- h+1 } nu<-c(1,2)[-mu]# for this direction we need to consider that the distance across the boundary# could be smaller h12<-min( (n[nu]-p[nu])^2, (1+n[nu]-p[nu])^2 )return(h^2+ h12)}forward_neigh_2d<-function(p,mu){ N<-dim(p)[1] res<-list() res$id<-array(dim =c( N,2 ))# id first neig, id second neigh res$h<-array(dim =c( N,2 ))# h first neig, h second neighfor( i inseq_along(p[,1])){ count <-1 neig<-list() neig$h<-rep(i, N-1) neig$id<-rep(i, N-1)for( j inseq_along(p[,1])[-i]){ neig$h[count]<-compute_distance_forward(p[i,],p[j,],mu) neig$id[count]<- j count<- count +1 }# print(neig$h)# print(neig$id) found<-0# print((neig$h))# print((neig$id))for (j inorder(neig$h)){# cat(j, found,"\n")if (found==0){# cat("found =0 \n") res$id[i,1]<-neig$id[j] res$h[i,1]<-neig$h[j] found <- found +1next }elseif (found==1){# cat("found =1 \n")#print(neig$h[j])#print(neig$id[j]) v1<-p[res$id[i,1], ]- p[i,] v2<-p[neig$id[j],]- p[i,] cth<-sum(v1*v2)/sqrt(sum(v1^2) *sum(v2^2))if( (1-abs(cth))>1e-9){ res$id[i,2]<-neig$id[j] res$h[i,2]<-neig$h[j] fount <- found +1break } }else{break } } }return(res)}nxp<-forward_neigh_2d(p,1)nxp
compute_distance2_wrapped<-function(p,n){# for the direction mu we compute the distance always forward # h1<-min( (n[1]-p[1])^2, (1+n[1]-p[1])^2 )# h2<-min( (n[2]-p[2])^2, (1+n[2]-p[2])^2 ) h1<- (n[1]-p[1])^2 h2<- (n[2]-p[2])^2return(h1+ h2)}wrap<-data.frame("x"=c(0,0,0,1,1,1,-1,-1,-1) , "y"=c(0,1,-1,0,1,-1,0,1,-1) )check_point_inside_wrapped<-function(center, r2 ,p){for (w inseq_along(wrap[,1])){ d2<-compute_distance2_wrapped( p+wrap[w,], center )if(d2<r2) return(TRUE) }return(FALSE)}#c(1,3,4), c(2,5,6)sort_links<-function(links, p1, p2 ){ to_switch<-which(links[,p1[1]]>links[,p2[1]])for (i inseq(p1)){ tmp <- links[to_switch, p2[i]] links[to_switch, p2[i]] <- links[to_switch, p1[i]] links[to_switch, p1[i]] <- tmp }return(links)}library(dplyr) #use to remove doublers
Attaching package: 'dplyr'
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
remove_doublers<-function(links){#df<- links[,c(1,2)]#df$l2<- (links$p1x-links$p2x)^2+(links$p1y-links$p2y)^2 o<-order(links$p1) links<- links[o,]#links<-distinct(links,p1,p2,l2,.keep_all = TRUE) rem<-c()for (i inc(1:(length(links$l2)-1) ) ){for (j inc((i+1):length(links$l2))){if (links$p1[i]==links$p1[j])if (links$p2[i]==links$p2[j])if (abs(links$l2[i]-links$l2[j])<1e-8){ rem<-c(rem,j) } } } links<-links[-rem,]#df<-distinct(df)return(links)}Delaunay_triangulation<-function(p){links<-data.frame("p1"=c(0), "p2"=c(0), "p1x"=c(0), "p1y"=c(0), "p2x"=c(0), "p2y"=c(0), "p3"=c(0), "p3x"=c(0), "p3y"=c(0) )#, "cx"=c(0), "cy"=c(0), "r"=c(0))links<-links[-1,]# we will compute the circle using the form# x^2+y^2+Ax+By+C=0# (x1 y1 1) (A) = -(x1^2+y1^2)# (x2 y2 1) (B) = -(x2^2+y2^2)# (x3 y3 1) (C) = -(x3^2+y3^2)# then we relare to # (x-cx)^2+(y-cy)^2=r^2M<-array(dim =c( 3,3 ))M[1,3]<-1M[2,3]<-1M[3,3]<-1b<-c(1,1,1)point_list<-seq_along(p[,1])point_list<-c(1: (length(p[,1])*9) )wp<- ppp<-ppp[,1]<-pp[,1]-1pp[,2]<-pp[,2]-1wp<-rbind(wp,pp )pp<-ppp[,2]<-pp[,2]-1wp<-rbind(wp,pp )pp<-ppp[,1]<-pp[,1]-1wp<-rbind(wp,pp )pp<- ppp[,2]<-pp[,2]+1wp<-rbind(wp,pp )pp<-ppp[,1]<-pp[,1]+1pp[,2]<-pp[,2]+1wp<-rbind(wp,pp )pp<-ppp[,1]<-pp[,1]+1wp<-rbind(wp,pp )pp<-ppp[,1]<-pp[,1]+1pp[,2]<-pp[,2]-1wp<-rbind(wp,pp )pp<-ppp[,1]<-pp[,1]-1pp[,2]<-pp[,2]+1wp<-rbind(wp,pp )wpm<-ppp<-ppp[,1]<-pp[,1]-1pp[,2]<-pp[,2]-1wpm<-rbind(wpm,pp )pp<-ppp[,2]<-pp[,2]-1wpm<-rbind(wpm,pp )pp<-ppp[,1]<-pp[,1]-1wpm<-rbind(wpm,pp )#we do not need to check the last one#point_list<-point_list[-length(point_list)]#for (i in point_list){for (i inseq_along(p[,1])){for (j inseq_along(wpm[,1])[-i] ){for (k inseq_along(wp[,1])[c(-i,-j)] ){ M[1,1]<-p[i,1] M[2,1]<-as.numeric(wp[j,1]) M[3,1]<-as.numeric(wp[k,1]) M[1,2]<-p[i,2] M[2,2]<-as.numeric(wp[j,2]) M[3,2]<-as.numeric(wp[k,2]) b[1]<--(M[1,1]^2+M[1,2]^2) b[2]<--(M[2,1]^2+M[2,2]^2) b[3]<--(M[3,1]^2+M[3,2]^2)if (abs(det(M)) <1e-6) next x<-solve(M,b) center<-c(-x[1]/2, -x[2]/2) r2<-center[1]^2+center[2]^2-x[3]# check if there are point inside the circle inside=FALSEfor (l inseq_along(wp[,1])[c(-i,-j,-k)]){# check all point since the wrapped can be inside d2<-compute_distance2_wrapped( wp[l,], center )if(d2<r2) inside<-TRUEif(inside) break }if(!inside){# we found a cluster of point jj<-((j-1)%%length(p[,1]))+1 kk<-((k-1)%%length(p[,1]))+1 links[length(links[,1])+1, c(1,2)]<-c(i,jj ) links[length(links[,1]), 3]<-M[1,1] links[length(links[,1]), 4]<-M[1,2] links[length(links[,1]), 5]<-M[2,1] links[length(links[,1]), 6]<-M[2,2] links[length(links[,1]), 7]<-kk links[length(links[,1]), 8]<-M[3,1] links[length(links[,1]), 9]<-M[3,2]# links[length(links[,1]), c(7,8,9)]<- c(center[1],center[2], sqrt(r2)) # ij<-TRUE# } links[length(links[,1])+1, c(1,2)]<-c(i,kk) links[length(links[,1]), 3]<-M[1,1] links[length(links[,1]), 4]<-M[1,2] links[length(links[,1]), 5]<-M[3,1] links[length(links[,1]), 6]<-M[3,2] links[length(links[,1]), 7]<-jj links[length(links[,1]), 8]<-M[2,1] links[length(links[,1]), 9]<-M[2,2]### links[length(links[,1]), c(7,8,9)]<- c(center[1],center[2], sqrt(r2))## no need to link j with k it will be found later# if(j<=length(p[,1]) || k<=length(p[,1]) ){# links[length(links[,1])+1, c(1,2)]<-c(jj,kk)# links[length(links[,1]), 3]<-M[2,1]# links[length(links[,1]), 4]<-M[2,2]# links[length(links[,1]), 5]<-M[3,1]# links[length(links[,1]), 6]<-M[3,2]# ### links[length(links[,1]), c(7,8,9)]<- c(center[1],center[2], sqrt(r2))# } } } }}links<-sort_links(links,c(1,3,4),c(2,5,6))links<-sort_links(links,c(1,3,4),c(7,8,9))links<-sort_links(links,c(2,5,6),c(7,8,9))links$l2<- (links$p1x-links$p2x)^2+(links$p1y-links$p2y)^2links$l13<- (links$p1x-links$p3x)^2+(links$p1y-links$p3y)^2links$l23<- (links$p2x-links$p3x)^2+(links$p2y-links$p3y)^2return(links)}
a faster version that rely on the fact that the noise on the point is less then the distance, thus adding noise to the point do not change the neighbor structure
We want to write a discrete version of the equation
\nabla^2 u=f
To do so we project both sides to a set of function h_i
\int (\nabla^2 u) h_i dV=\int f h_i dV
\tag{1}
since our lattice is triangulated we can define write our integral as a sum over the integral on the triangles \Delta_i
\int f dV=\sum_{i\in \Delta}\int_{\Delta_i} f dV
we can write f=\sum f_i h_i , if h_i where a basis then it will be exact but in our case it is an approximation. The h_i are chosen to be 1 in the point i and zero in the all the other point. They go linear from 0 to 1 thus
h_i=ax+by+c
more specific given a trinagle with vertices i,j,k the function h_i satify \begin{align}
h_i(i)&=ax_i+by_i+c=1\\
h_i(j)&=ax_j+by_j+c=0\\
h_i(k)&=ax_k+by_k+c=0\,.
\end{align} To approximate the right had side of (Equation 1) we tried:
writing f=\sum_i f_i h_i\begin{equation}
\int f h_i dV=\sum_{ \Delta_i } \sum_{j=0}^N f_j\int_{\Delta_i} h_ih_jdV =M_{ij}f_j
\end{equation} where \sum_{\Delta_i} is the summ of all triangles having i as vertex and the integral can be approximate with the midpoint rule for triangles \begin{equation}
\int_{\Delta_i} h_ih_jdV =\frac{Area}{3}\left[h_ih_j\big|_{(\frac{\bar x_i+\bar x_j}{2})}+h_ih_j\big|_{(\frac{\bar x_i+\bar x_k}{2})}+h_ih_j\big|_{(\frac{\bar x_j+\bar x_k}{2})} \right]
\end{equation}
Important
This works with regular Delanuy but not in the generic case
approximate f\approx f_i inside the triangle, then the integral can be done exactly and is equal to the volume of the pyramid \begin{equation}
\int f h_i dV = f_i \frac{1}{3}\sum_{ \Delta_i } Area_{ \Delta_i }
\end{equation} in this case M is diagonal.
Important
This works with regular Delanuy but still does not in the generic case, to make it work with a generic triangularization we found that we need to replace M with its average value M=\bar{M}1\!\!1
To approximate the left hand side of (Equation 1) we write u=\sum_j u_jh_j integrate by part \begin{align}
\int (\nabla^2 u) h_i dV&=\sum_{i \in \Delta} \sum_{j=0}^N u_j\int_{\Delta_i}(\nabla^2h_j)h_idV
\nonumber\\
&=\sum_{i \in \Delta} \sum_{j=0}^N u_j\int_{\Delta_i}-(\nabla h_j)(\nabla h_i)dV+
\sum_{i \in \Delta} \sum_{j=0}^N u_j\int_{\Delta_i}\nabla[(\nabla h_j) h_i]dV
\nonumber\\
&=\sum_{i \in \Delta} \sum_{j=0}^N u_j\int_{\Delta_i}-(\nabla h_j)(\nabla h_i)dV+
\sum_{i \in \Delta} \sum_{j=0}^N u_j\int_{S_{\Delta_i}}[(\nabla h_j) h_i]\cdot n \,dS
\nonumber\\
&=\sum_{i \in \Delta} \sum_{j=0}^N u_j\int_{\Delta_i}-(\nabla h_j)(\nabla h_i)dV
=L_{ij}u_j
\end{align} from the second to the third line we use the divergence theorem to tranform the integral over the volume as the flux on the surface, since our lattice is triangulated the flux on the surface cancel between the neighbor triangles.
Thus our discretized version of the Laplace equation become
L_{ij}u_j=M_{ij}f_j
and we can write a discretized version of the Laplatian as
\nabla^2\to M^{-1}L
library(ggplot2)library(plotly)library(matlib)library(Matrix)compute_hi<-function(p,M,i){ b<-c(0,0,0) b[i]<-1 a<-solve(M,b) A<-det(M)/2return(a*sqrt(abs(A)) )}nablah_dot_nablah<-function(h1,h2){return(h1[1]*h2[1]+h1[2]*h2[2] )}func_h<-function(h, p ){return(h[1]*p[1]+h[2]*p[2]+h[3])}midpoint_triangle<-function(hi, hj, v1,v2,v3){return(func_h(hi, (v1+v2)/2 )*func_h(hj, (v1+v2)/2 )+func_h(hi, (v1+v3)/2 )*func_h(hj, (v1+v3)/2 )+func_h(hi, (v2+v3)/2 )*func_h(hj, (v2+v3)/2 ) )}## type is an integer that defines how the rhs of laplace equation is computed## 0 trapezoidal rule## 1 trapezoidal rule + average of the diagonal of M## 2 gauss quadrature 3 point compute_laplacian<-function(p,triangles, type=1){#print(triangles) L<-array(rep(0,length(p[,1])^2), dim =c( length(p[,1]),length(p[,1]) ))#Matrix(nrow = length(p[,1]), ncol = length(p[,1]), data = 0, sparse = TRUE)# M<-array(rep(0,length(p[,1])^2), dim =c( length(p[,1]),length(p[,1]) ))#Matrix(nrow = length(p[,1]), ncol = length(p[,1]), data = 0, sparse = TRUE)# for (i inseq_along(triangles[,1]) ){ v1<-c(triangles[i,4], triangles[i,5]) v2<-c(triangles[i,6], triangles[i,7]) v3<-c(triangles[i,8], triangles[i,9]) Tri<-t(array( c( v1[1],v1[2],1, v2[1],v2[2],1, v3[1],v3[2],1),dim =c( 3,3 )) ) Area <-abs(det(Tri))/2 h<-array( dim =c( 3,3 ))# cat("triangle out: ", triangles[i,][[1]], triangles[i,][[2]], triangles[i,][[3]],"\n" ) h[1,] <-solve(Tri, c(1,0,0))# compute_hi(p, Tri,1) h[2,] <-solve(Tri, c(0,1,0))#compute_hi(p, Tri,2) h[3,] <-solve(Tri, c(0,0,1))#compute_hi(p, Tri,3)for (j inc(1,2,3)){if (type==2){for (k inc(1,2,3)){ M[triangles[i,j],triangles[i,k] ]<- M[triangles[i,j],triangles[i,k]] + (Area/3)*midpoint_triangle(h[j,], h[k,], v1, v2, v3 ) } }else{ M[triangles[i,j],triangles[i,j] ] <- M[triangles[i,j],triangles[i,j] ] + (Area/3) } }for (j inc(1,2,3)){for (k inc(1,2,3)){ L[triangles[i,j], triangles[i,k] ] <- L[triangles[i,j], triangles[i,k] ]-nablah_dot_nablah(h[j,], h[k,])*Area } } }if(type==2|| type ==0 ){ L <--inv(M) %*% L #we add the -1 since we want (-i\nabla)^2 }elseif (type==1){ L<-- (1/mean(diag(M)))*L }return(L)}#sqN<- 5#r<- 1e-6#p <- noise_on_the_regular_lattice_2d(sqN,r/sqN)# f <- sin(2*pi*p[,1])+cos(2*pi*p[,2])# links <- Delaunay_triangulation_fast(p)# links1 <- find_triangles(links)# L2 <- compute_laplacian(p,links1)# plot_Delaunay_triangulation(p,links1)
Test application on a function
Consider the function
f(x,y)=\sin (2\pi x)+\cos (2\pi y)
with the lapacian
\nabla^2f(x,y)=-(2\pi)^2[\sin (2\pi x)+\cos (2\pi y)]=-(2\pi)^2f(x,y)
We define the deviation between exact derivative and the discrete as (we want to implemente (-i\nabla)^2=-\nabla^2:
d=\frac{1}{N}\sum_i^N |\nabla^2_{latt}f-\nabla^2_{cont}f |
On a regular lattice of size 5x5 we have |\nabla^2_{latt}f-\nabla^2_{cont}f|\approx\epsilon\approx\frac{1}{\sqrt{N}} thus d\sim \frac{1}{\sqrt{N}}
On a regular lattice there is an ambiguity in defining the Delaunay triangularization since we have 4 point on the circumference. We solve the ambiguity definin it in this way
Ns<-c(8,10,12)dfs<-data.frame("id"=c(0), "l"=c(0), "N"=c(0), "type"=c(0), "r"=c(0))dfs<-dfs[-1,]r<-1e-6for (N in Ns){ p <-noise_on_the_regular_lattice_2d(N,r/N) N<-N*N links <-Delaunay_triangulation_fast(p) links1<-find_triangles(links) L2<-compute_laplacian(p,links1) l <-eigen(L2)$values[(N):1] # divide by N to remove the lattice spacing dfs<-rbind(dfs,data.frame("id"=c(1:N)-1, "l"=l, "N"=rep(N,N),"type"=rep("L2",N) , "r"=rep(r,N)) )}
Ns<-c(8,10,12,14)dfrs<-data.frame("id"=c(0), "l"=c(0), "N"=c(0), "type"=c(0), "r"=c(0))dfrs<-dfrs[-1,]r<-1e-6for (N in Ns){ p <-noise_on_the_regular_lattice_2d(N,r/N) N<-N*N links <-Delaunay_triangulation_regular(p) L2<-compute_laplacian(p,links) l <-eigen(L2)$values[(N):1] # divide by N to remove the lattice spacing dfrs<-rbind(dfrs,data.frame("id"=c(1:N)-1, "l"=l, "N"=rep(N,N),"type"=rep("L2",N) , "r"=rep(r,N)) )}