Improved space-time plot. Fixed Domain constructor for sf objects.

fdaPDE · Nov 28, 2023 · 73633ca · 73633ca
1 parent 741dc2b
commit 73633ca
Show file tree

Hide file tree

Showing 7 changed files with 51 additions and 33 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -1,9 +1,13 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method("st_crs<-",DomainObject)
 S3method(plot,FunctionObject)
 S3method(plot,MeshObject)
+S3method(st_as_sf,DomainObject)
+S3method(st_as_sfc,DomainObject)
+S3method(st_as_sfc,MeshObject)
+S3method(st_crs,DomainObject)
 export("%X%")
-export("st_crs<-.DomainObject")
 export(Domain)
 export(Function)
 export(FunctionSpace)
@@ -15,10 +19,6 @@ export(dot)
 export(grad)
 export(laplace)
 export(read_mesh)
-export(st_as_sf.DomainObject)
-export(st_as_sfc.DomainObject)
-export(st_as_sfc.MeshObject)
-export(st_crs.DomainObject)
 exportMethods("*")
 exportMethods("+")
 exportMethods("-")

diff --git a/R/domain.R b/R/domain.R
@@ -176,7 +176,13 @@ setMethod("Domain", signature="sfc", function(x){
         node_group_i <-matrix(0, nrow(nodes[id_edge_sub,]),ncol=1)
         for(e in 1:nrow(edge_sub_i)){
           if( nodes_ring[edge_sub_i[e,1],ncol(coords)] & nodes_ring[edge_sub_i[e,2],ncol(coords)] ){ # sono entrambi sul bordo
-            edge_bd_i[e] <- 1
+            if(!st_within( st_centroid(st_linestring(as.matrix(rbind(nodes_ring[edge_sub_i[e,1],1:2], # controllo che il punto medio 
+                                                         nodes_ring[edge_sub_i[e,2],1:2])))),         # del lato NON giaccia all'interno -> nodi di bordo che collegano due lati
+                          st_union(st_geometry(x)),                                                   # NB. sembra rallentare
+                          sparse=F)){
+              edge_bd_i[e] <- 1 
+
+            } 
             edge_group_i[e] <- sub_id
             #edge_BC_i[e] <- "neumann"
             node_group_i[edge_sub_i[e,1]] <- e

diff --git a/R/function_space.R b/R/function_space.R
@@ -122,7 +122,7 @@ plot.FunctionObject <- function(x, ...){
   K=x$FunctionSpace$mesh$get_elements()[,3]-1
   fig<- plot_ly(plot_data, x=~X, y=~Y, z=~coeff,
           i = I, j = J, k = K,
-          intensity=~coeff,color = ~coeff, type="mesh3d", 
+          intensity=~coeff, color=~coeff, type="mesh3d", 
           colorbar=list(title=""), ...) %>%
       layout(scene = list(
         aspectmode = "data", 
@@ -138,15 +138,16 @@ plot.FunctionObject <- function(x, ...){
   }else{
     plot_data <- data.frame(X=rep(x$FunctionSpace$mesh$get_nodes()[,1], times=length(times)), 
                             Y=rep(x$FunctionSpace$mesh$get_nodes()[,2], times=length(times)),
+                            Z=rep(0, times=length(times)),
                             coeff=as.vector(x$coeff[1:nrow(x$FunctionSpace$mesh$get_nodes()),]),
-                            times = rep(times, each=nrow(x$FunctionSpace$mesh$get_nodes())))
+                            times = round(rep(times, each=nrow(x$FunctionSpace$mesh$get_nodes())),3))
     limits = c(min(x$coeff), max(x$coeff))
     I=x$FunctionSpace$mesh$get_elements()[,1]-1
     J=x$FunctionSpace$mesh$get_elements()[,2]-1
     K=x$FunctionSpace$mesh$get_elements()[,3]-1
-    fig<- plot_ly(plot_data, x=~X, y=~Y, z=~coeff, frame=~times,
+    fig<- plot_ly(plot_data, x=~X, y=~Y, z=~Z, frame=~times,
                   i = I, j = J, k = K, cmin = limits[1], cmax=limits[2],
-                  intensity=~coeff,color = ~coeff, type="mesh3d",
+                  intensity=~coeff, color=~coeff, type="mesh3d",
                   colorbar=list(title=""), ...) %>%
         layout(scene = list(
           aspectmode = "data", 
@@ -155,12 +156,12 @@ plot.FunctionObject <- function(x, ...){
           yaxis = list(
             title = '', showgrid = F, zeroline = F, showticklabels = F),
           zaxis = list(
-            title = '', showgrid = F, zeroline = F, showticklabels = F)),
+            title = '', showgrid = F, zeroline = F, showticklabels = F),
           camera = list(
-            eye = list(x = 1.25, y = -1.25, z = 1.25))) %>%
+            eye = list(x = 0, y = -0.01,  z = 1.25))), dragmode="zoom") %>%
         colorbar(len = 1, title="") %>%
-        animation_opts(frame=5) %>% 
         animation_slider(currentvalue = list(prefix ="t = "))
+        #animation_opts(frame=5) %>% 
   }
   fig 
 }
@@ -220,8 +221,8 @@ setMethod("contour", signature=c(x="FunctionObject"), function(x, ...){
                    colorbar=list(title=""), ...) %>%
       layout(xaxis = list(title = "", showgrid=F, zeroline=F, ticks="", showticklabels=F),
              yaxis = list(title = "", showgrid=F, zeroline=F, ticks="", showticklabels=F)) %>%
-      animation_opts(frame=5) %>% 
-      animation_slider(currentvalue = list(prefix ="t = "))
+      animation_slider(currentvalue = list(prefix ="t = ")) #animation_opts(frame=5) %>% 
+
   }
   fig
 }

diff --git a/tests/build_domain.R b/tests/build_domain.R
@@ -121,7 +121,24 @@ st_as_sf(domain)
 tmp <- st_as_sf(domain)
 mesh <- build_mesh(domain, maximum_area = 0.0025, minimum_angle = 20)
 plot(st_as_sfc(mesh))
-points(mesh$get_nodes()[mesh$get_boundary()==1,], pch=16, col="red") # da sistemare ?
+points(mesh$get_nodes()[mesh$get_boundary()==1,], pch=16, col="red") 
+
+st_intersection(st_centroid(st_linestring(nodes[c(5,6),])), geom_sfc)
+
+st_contains(st_geometry(geom_sfc), st_centroid(st_linestring(nodes[c(5,6),])),sparse = F)
+
+st_within(st_centroid(st_linestring(nodes[c(5,6),])), 
+          st_union(st_geometry(geom_sfc)), 
+          sparse=F)
+
+st_within(st_centroid(st_linestring(nodes[c(1,2),])), 
+          st_union(st_geometry(geom_sfc)), 
+          sparse=F)
+
+plot(st_intersection(geom_sfc, st_linestring(nodes[c(5,6),])), col="red")
+plot(st_difference(geom_sfc, st_linestring(nodes[c(5,6),])), col="red")
+st_intersection(st_difference(geom_sfc, st_linestring(nodes[c(5,6),])),
+                st_linestring(nodes[c(5,6),]))
 # ------------------------------------------------------------------------------
 
 # reading pslg object ----------------------------------------------------------

diff --git a/tests/parabolic.R b/tests/parabolic.R
@@ -80,19 +80,13 @@ dim( u$eval_at(points) )
 max(abs( u$eval_at(points) - as.matrix(exact_solution(points))))
 # 
 # ## plot solution 
-# options(warn=-1)
-plot(u) %>% hide_colorbar()
+options(warn=-1)
+plot(u) 
 
-# scene = list(camera = list(eye = list(x = 0, y = 0, z = 1.25)),
-#              xaxis = list(autorange="reversed"),
-#              yaxis = list(autorange="reversed"))
-# 
-# plot_tmp <- plot(u) %>% 
-#   layout(dragmode=FALSE, scene=scene,aspectmode="cube") %>%
-#   hide_colorbar() %>% 
-#   config(displayModeBar =FALSE)
-# 
-# plot_tmp
+plot(u) %>% layout(scene =list(camera=list(eye=list(z=1))))
+plot(u, showscale=FALSE) # no color bar :)
+
+plot(u)
 
 # contour
 contour(u)
diff --git a/vignettes/Introduction.Rmd b/vignettes/Introduction.Rmd
@@ -56,8 +56,8 @@ p <- pslg(P=rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 1)),
           S=rbind(c(1, 2), c(2, 3), c(3, 4), c(4,1)))
 unit_square <- triangulate(p, a = 0.00125, q=30)
 domain <- Mesh(unit_square)
-xrange <- range(domain$nodes()[,1])
-yrange <- range(domain$nodes()[,2])
+xrange <- range(domain$get_nodes()[,1])
+yrange <- range(domain$get_nodes()[,2])
 Nx <- 40
 Ny <- 40
 eval_x <- seq(from=xrange[1], to=xrange[2], length.out=Nx)
@@ -193,7 +193,7 @@ mesh_aux <- Mesh(tpp)
 Vh_aux <- FunctionSpace(mesh_aux)
 hat_function <- Function(Vh_aux)
 
-coeff_aux <- rep(0, times=nrow(mesh$nodes()))
+coeff_aux <- rep(0, times=nrow(mesh$get_nodes()))
 coeff_aux[5] = 1
 
 hat_function$set_coeff(as.matrix(coeff_aux))

diff --git a/vignettes/Parabolic.Rmd b/vignettes/Parabolic.Rmd
@@ -33,7 +33,7 @@ $$
 
 
 where $u_0 = \sin( 2\pi x)\sin(2\pi y)$ is the initial condition, $f(x,y) = 8\pi^2 \sin( 2\pi x)\sin(2\pi y) e^{-t}$ is the forcing term and $\partial \Omega$ is the boundary of $\Omega$ where we have prescribed homogeneous Dirichelet boundary condition for all $t$ grater than $0$.
-The exact solution of the previous problem is $u_{ex}(x,y) = \sin(2\pi x) \sin(2 \pi y) e^{-t}$ whose trend is shown in the following.
+The exact solution of the previous problem is $u_{ex}(x,y) = \sin(2\pi x) \sin(2 \pi y) e^{-t}$. The following window wraps all the steps needed to solve the problem.
 
 <!-- ```{r, echo=FALSE} -->
 <!-- exact_solution <- function(points,t){ -->
@@ -126,5 +126,5 @@ pde$set_dirichletBC(dirichletBC)
 pde$solve()
 
 ## plot solution
-contour(u)
+plot(u)
 ```