diff --git a/NAMESPACE b/NAMESPACE
index 790bcc3..456685d 100644
--- 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("-")
@@ -42,4 +42,5 @@ importFrom(sf,st_linestring)
 importFrom(sf,st_point)
 importFrom(sf,st_polygon)
 importFrom(sf,st_sfc)
+importFrom(sf,st_within)
 useDynLib(femR)
diff --git a/R/domain.R b/R/domain.R
index e847f98..5d0dd53 100644
--- a/R/domain.R
+++ b/R/domain.R
@@ -5,17 +5,6 @@
                             time_interval = "numeric",
                             crs ="ANY"
                            )
-                           #,
-                           # methods=c(                                 # on -> list(sub_id=sub_id, edges_id=) local edge id(s)   
-                           #    set_boundary_type = function(on, type){ # "neumann", "dirichlet"
-                           #      if(is.null(on$sub_id)) sub_id <- 1
-                           #      edges_id <- on$edges_id 
-                           #      geometry[[sub_id]]$BC[edges_id] <<- rep(type, times=length(edges_id))
-                           #      nodes_id <- unique(as.vector(geometry[[sub_id]]$edges[edges_id,]))
-                           #      #geometry[[sub_id]]$nodes_group[nodes_id] <<- rep(type, length(nodes_id))
-                           #    }
-                           # )
-                          
 )
 
 #' S4 class representing a spatial (spatio-temporal) domain
@@ -28,6 +17,7 @@
 #' @rdname DomainObject
 #' @importFrom sf NA_crs_
 #' @importFrom sf st_crs
+#' @importFrom sf st_within
 #' @export 
 setGeneric("Domain", function(x) standardGeneric("Domain"))
 
@@ -38,10 +28,6 @@ setMethod("Domain", signature = "list", function(x){
     if(is.null(x$nodes_boundary)) x$nodes_boundary <- rep(1, times=nrow(x$nodes))
     if(is.null(x$edges_boundary)) x$edges_boundary <- rep(1, times=nrow(x$edges))
     if(is.null(x$edges_ring))     x$edges_ring <- rep(1, times = nrow(x$edges))
-    # if(is.null(x$BC)){
-    #   x$BC <- rep(NA, length=length(x$edges_boundary))
-    #   x$BC[which(x$edges_boundary==1)] = "neumann"
-    # } 
     
     geometry <- vector(mode="list", length=1L)
     geometry[[1]] <- list(nodes=x$nodes, nodes_group=x$nodes_group, nodes_boundary=x$nodes_boundary,
@@ -50,7 +36,6 @@ setMethod("Domain", signature = "list", function(x){
     
     coords <- data.frame(x=x$nodes[,1], y=x$nodes[,2], boundary=x$nodes_boundary)
   crs = NA
-  #if(!is.null(x$crs)) crs <- st_crs(x$crs)
   .DomainCtr(geometry = geometry, time_interval = vector(mode="numeric", length = 0), 
              coords=coords, crs = crs)
 })
@@ -68,8 +53,6 @@ setMethod("Domain", signature = "pslg", function(x){
   
   nodes_group <- rep(1, times=nrow(nodes))
   edges_group <- rep(1, times=nrow(edges))
-  #BC <- rep(NA, length=length(edges_boundary))
-  #BC[which(edges_boundary==1)] = "neumann"
   edges_ring <- rep(1, times=nrow(edges))
   
   geometry <- vector(mode="list", length=1L)
@@ -86,59 +69,54 @@ setMethod("Domain", signature = "pslg", function(x){
 #' @rdname DomainObject
 setMethod("Domain", signature="sfc", function(x){
   if(length(x) == 1) x <- st_cast(x, to="MULTIPOLYGON") 
-  coords <- st_coordinates(x)
-  #coords <- coords[-nrow(coords),]
-  coords <- cbind(coords, 
-                  rep(0, times=nrow(coords)), # id 
-                  rep(0, times=nrow(coords))) # boundary (1,0)
-  coords <- as.data.frame(coords)
-  colnames(coords)[(ncol(coords)-1):ncol(coords)] <- c("id", "boundary")
-  storage.mode(coords$id) <- "integer"; storage.mode(coords$boundary) <- "integer" 
-  pts_list <- vector(mode="list", length=nrow(coords))
-  for(i in 1:nrow(coords)){
-    pts_list[[i]] <- st_point(as.numeric(coords[i,1:2]))
-  }
+  st_coords <- st_coordinates(x)
+  
+  st_coords <- cbind(st_coords, 
+                  rep(0, times=nrow(st_coords)), # id 
+                  rep(0, times=nrow(st_coords))) # boundary (1,0)
+  st_coords <- as.data.frame(st_coords)
+  colnames(st_coords)[(ncol(st_coords)-1):ncol(st_coords)] <- c("id", "boundary")
+  storage.mode(st_coords$id) <- "integer"; storage.mode(st_coords$boundary) <- "integer" 
+  
+  # nodes to st_point ! (FAST)
+  pts_list <- lapply(as.list(as.data.frame(t(st_coords[,1:2]))), st_point)
   pts_sfc <- st_sfc(pts_list)
   
   # setto indici dei nodi al bordo
   coords_bd <- st_coordinates(st_cast(st_union(st_geometry(x)),
                                       to="LINESTRING"))
-  pts_bd_list <- vector(mode="list", length=nrow(coords_bd))
-  for(i in 1:nrow(coords_bd)){
-    pts_bd_list[[i]] <- st_point(coords_bd[i,1:2])
-  }
+  
+  pts_bd_list <- lapply(as.list(as.data.frame(t(coords_bd[,1:2]))), st_point)
   pts_bd_sfc <- st_sfc(pts_bd_list)
   dist_bd <- st_distance(pts_sfc, pts_bd_sfc)
   
   num_nodes <- 1
   for(i in 1:nrow(coords_bd)){
     id_bd <- which(dist_bd[,i] < 100 *.Machine$double.eps )
-    if(any( coords[id_bd,(ncol(coords)-1)] != 0)){
-      global_id <- which(coords[id_bd,(ncol(coords)-1)] != 0)[1]
-      coords[id_bd,(ncol(coords)-1)] <-  rep(global_id, times=length(id_bd))
+    if(any( st_coords[id_bd,(ncol(st_coords)-1)] != 0)){
+      global_id <- which(st_coords[id_bd,(ncol(st_coords)-1)] != 0)[1]
+      st_coords[id_bd,(ncol(st_coords)-1)] <-  rep(global_id, times=length(id_bd))
     }else{
-      coords[id_bd,(ncol(coords)-1)] <- rep(num_nodes, times=length(id_bd))
+      st_coords[id_bd,(ncol(st_coords)-1)] <- rep(num_nodes, times=length(id_bd))
       num_nodes <- num_nodes + 1  
     }
-    coords[id_bd, ncol(coords)] <- rep(1, times=length(id_bd))
+    st_coords[id_bd, ncol(st_coords)] <- rep(1, times=length(id_bd))
   }
   
   for(sub_id in 1:length(st_geometry(x))){
     coords_sub_id <- st_coordinates(st_geometry(x)[sub_id])
-    pts_sub_id_list <- vector(mode="list", length=nrow(coords_sub_id))
-    for(j in 1:nrow(coords_sub_id)){
-      pts_sub_id_list[[j]] <- st_point(coords_sub_id[j,1:2])
-    }
+    pts_sub_id_list <- lapply(as.list(as.data.frame(t(coords_sub_id[,1:2]))), st_point)
+    
     pts_sub_id_sfc <- st_sfc(pts_sub_id_list)
     dist_matrix_sub_id <- st_distance(pts_sfc, pts_sub_id_sfc)
     for(j in 1:nrow(coords_sub_id)){
       id <- which(dist_matrix_sub_id[,j] < 100* .Machine$double.eps)
-      if(all(coords[id,(ncol(coords)-1)] != 0 )){ next }
-      if( any( coords[id,(ncol(coords)-1)] != 0) ){
-        global_id <- which(coords[id,(ncol(coords)-1)] != 0)[1]
-        coords[id,(ncol(coords)-1)] <-  rep(global_id, times=length(id))
+      if(all(st_coords[id,(ncol(st_coords)-1)] != 0 )){ next }
+      if( any( st_coords[id,(ncol(st_coords)-1)] != 0) ){
+        global_id <- which(st_coords[id,(ncol(st_coords)-1)] != 0)[1]
+        st_coords[id,(ncol(st_coords)-1)] <-  rep(global_id, times=length(id))
       }else{
-        coords[id,(ncol(coords)-1)] <- rep(num_nodes, times = length(id))
+        st_coords[id,(ncol(st_coords)-1)] <- rep(num_nodes, times = length(id))
         num_nodes = num_nodes + 1
       }
     }
@@ -147,79 +125,83 @@ setMethod("Domain", signature="sfc", function(x){
   geometry <- vector(mode="list", length =length(st_geometry(x)))
   
   for(sub_id in 1:length(st_geometry(x))){
-    idx <- which(coords[,5]==sub_id)  
-    nodes <- unique(coords[idx,])
-    n_main_ring <- diff(range(nodes[,4])) + 1 # number of main ring!!! see st_coordinates
-    edge <- matrix(nrow=0, ncol=2)
-    edge_sub <- matrix(nrow=0, ncol=2)
-    edge_group <- matrix(NA, nrow=0, ncol=1)
-    edge_bd <- matrix(nrow=0, ncol=1)
-    #edge_BC <- matrix(NA,nrow=0, ncol=1)
-    edge_ring <- matrix(nrow=0, ncol=1)
-    node_group <- matrix(nrow=0, ncol=1)
-    for(ring in 1:n_main_ring){
-      idx_ring <- which(nodes[,4] == ring)
-      nodes_ring <- nodes[idx_ring,]
+      
+    nodes <- unique(st_coords[which(st_coords[,5]==sub_id),])
+    n_rings <- diff(range(nodes[,4])) + 1 # number of "rings"! see st_coordinates,
+                                          # 1 main ring (external boundary), 
+                                          #>1 interior rings (holes!)
+    
+    edges <- matrix(nrow=0, ncol=2)        
+    
+    edges_group <- matrix(NA, nrow=0, ncol=1)
+    edges_bd <- matrix(nrow=0, ncol=1)
+    edges_ring <- matrix(nrow=0, ncol=1)
+    nodes_group <- matrix(nrow=0, ncol=1)
+    for(ring in 1:n_rings){
+
+      nodes_ring <- nodes[which(nodes[,4] == ring),]
       # in the current main ring
-      n_physical_lines <- diff(range(nodes_ring[,3])) + 1 
-      for(i in 1:n_physical_lines){
+      n_lines <- diff(range(nodes_ring[,3])) + 1 
+      for(i in 1:n_lines){
         id_edge_sub <- which(nodes_ring[,3] == i)
         # local
-        edge_sub_i <-  cbind(id_edge_sub,
+        edges_sub_i <-  cbind(id_edge_sub,
                              c(id_edge_sub[2:length(id_edge_sub)],id_edge_sub[1]))
         # global 
-        edge_i <- cbind(nodes_ring[edge_sub_i[,1],(ncol(coords)-1)], nodes_ring[edge_sub_i[,2],(ncol(coords)-1)])           
+        edges_i <- cbind(nodes_ring[edges_sub_i[,1],(ncol(st_coords)-1)], nodes_ring[edges_sub_i[,2],(ncol(st_coords)-1)])           
+        
+        edges_bd_i <- matrix(0, nrow=nrow(edges_sub_i),ncol=1)
+        edges_group_i <- matrix(NA, nrow=nrow(edges_sub_i), ncol=1)
+        
+        nodes_group_i <-matrix(0, nrow(nodes[id_edge_sub,]),ncol=1)
+        
+        centroids_list <- lapply( as.list(as.data.frame(t(edges_sub_i))), FUN=function(edge){
+          st_centroid(st_linestring(as.matrix(nodes_ring[edge,1:2])))
+        })
         
-        edge_bd_i <- matrix(0, nrow=nrow(edge_sub_i),ncol=1)
-        edge_group_i <- matrix(NA, nrow=nrow(edge_sub_i), ncol=1)
-        #edge_BC_i <- matrix(NA, nrow=nrow(edge_sub_i),ncol=1)
-        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
-            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
-            node_group_i[edge_sub_i[e,2]] <- e
+        is_inside <- st_within(st_sfc(centroids_list, crs=st_crs(x)$input), 
+                                    st_union(st_geometry(x)), sparse = F)
+        
+        for(e in 1:nrow(edges_sub_i)){
+          if( nodes_ring[edges_sub_i[e,1],ncol(st_coords)] & nodes_ring[edges_sub_i[e,2],ncol(st_coords)] ){ 
+            if(!is_inside[e]) edges_bd_i[e] <- 1 
+            
+            edges_group_i[e] <- sub_id
+            nodes_group_i[edges_sub_i[e,1]] <- e
+            nodes_group_i[edges_sub_i[e,2]] <- e
           } 
         }
-        if(i == 1){
-          edge_ring_i <- as.matrix(rep(ring, times=nrow(edge_sub_i)))
-        }else{
-          edge_ring_i <- as.matrix(rep((-i+1), times=nrow(edge_sub_i)))
+        
+        if(i == 1){ 
+          edge_ring_i <- as.matrix(rep(ring, times=nrow(edges_sub_i)))
+        }else{  
+          edge_ring_i <- as.matrix(rep((-i+1), times=nrow(edges_sub_i)))
         }
-        edge <- rbind(edge, edge_i)
-        edge_ring <- rbind(edge_ring, edge_ring_i)
-        edge_bd <- rbind(edge_bd, edge_bd_i)
-        edge_group <- rbind(edge_group, edge_group_i)
-        #edge_BC <- rbind(edge_BC, edge_BC_i)
-        edge_sub <- rbind(edge_sub, edge_sub_i)
-        node_group <- rbind(node_group, node_group_i)
+        edges <- rbind(edges, edges_i)
+        edges_ring <- rbind(edges_ring, edge_ring_i)
+        edges_bd <- rbind(edges_bd, edges_bd_i)
+        edges_group <- rbind(edges_group, edges_group_i)
+        #edge_sub <- rbind(edge_sub, edges_sub_i)
+        nodes_group <- rbind(nodes_group, nodes_group_i)
       }
     }
-    storage.mode(edge) <-"integer"
-    storage.mode(edge_bd) <-"integer"
-    storage.mode(nodes[,ncol(coords)]) <- "integer" 
-    geometry[[sub_id]] <- list(nodes = as.matrix(nodes[,1:2]), nodes_group = node_group, nodes_boundary = nodes[,ncol(coords)],
-                               edges = edge, edges_group = edge_group, edges_boundary = edge_bd,
-                               edges_ring = edge_ring)
+    storage.mode(edges) <-"integer"
+    storage.mode(edges_bd) <-"integer"
+    storage.mode(nodes[,ncol(st_coords)]) <- "integer" 
+    geometry[[sub_id]] <- list(nodes = as.matrix(nodes[,1:2]), nodes_group = nodes_group, nodes_boundary = nodes[,ncol(st_coords)],
+                               edges = edges, edges_group = edges_group, edges_boundary = edges_bd,
+                               edges_ring = edges_ring)
     
   }
   # da salvare
-  NODES <- unique(coords[order(coords[,(ncol(coords)-1)]),][,c(1:2,ncol(coords))])
-  NODES <- as.data.frame(NODES)
-  storage.mode(NODES$boundary) <- "integer"
+  coords <- unique(st_coords[order(st_coords[,(ncol(st_coords)-1)]),][,c(1:2,ncol(st_coords))])
+  coords <- as.data.frame(coords)
+  storage.mode(coords$boundary) <- "integer"
   
   crs <- NA
   if(!is.na(st_crs(x))) crs <- st_crs(x)$input
       .DomainCtr(geometry = geometry, time_interval = vector(mode="numeric", length = 0), 
-                coords=NODES, crs = crs)
+                coords=coords, crs = crs)
 })
 
 # setMethod("Domain", signature = "sfc", function(x){
@@ -261,7 +243,7 @@ setMethod("build_mesh", signature=c("DomainObject", "numeric", "numeric"),
     groups_id <- unique(domain$geometry[[sub_id]]$edges_ring)
     holes_id <- groups_id[ groups_id < 0 ]
     num_holes <- length(holes_id)
-    #holes <- matrix(0, nrow=num_holes, ncol=2)
+    
     if(num_holes > 0){
       for(i in 1:num_holes){
         edges_id <- which(domain$geometry[[sub_id]]$edges_ring==holes_id[i])
@@ -286,55 +268,12 @@ setMethod("build_mesh", signature=c("DomainObject", "numeric", "numeric"),
                S = S, SB = SB, H = H) 
   triangulation <- triangulate(p = pslg, a = maximum_area, q = minimum_angle)
   
-  # num_bd_nodes <- sum(triangulation$PB==1)
-  # bd_id <- which(triangulation$PB==1)
-  # pts_list <- list()
-  # for(i in 1:num_bd_nodes){
-  #   pts_list[[i]] <- st_point( t(as.matrix(triangulation$P[bd_id[i],1:2])))
-  # }
-  # 
-  # triangulation$PB <- matrix(0, nrow=nrow(triangulation$P),ncol=1)
-  # for(e in 1:nrow(S)){
-  #   if(SB[e] == 1 & BC[e] == "dirichlet"){
-  #     edge <- st_linestring(as.matrix(triangulation$P[S[e,],1:2]))
-  #     dirichlet_nodes <- st_intersects(edge, 
-  #                                      st_sfc(pts_list),sparse = FALSE)
-  #     dirichlet_id <- which(dirichlet_nodes[1,] == T)
-  #     triangulation$PB[ bd_id[dirichlet_id] ] = 1
-  #   }
-  # }
   res <- Mesh(triangulation)
   res$geometry <- domain$geometry
   res$time_interval <- domain$time_interval
   res$crs <- domain$crs
   return(res)
 })
-# setMethod("build_mesh", signature=c("DomainObject","numeric", "numeric"), 
-#           function(domain, maximum_area, minimum_angle){
-#             groups_id <- unique(domain$geometry$edges_ring)
-#             holes_id <- groups_id[ groups_id < 0 ]
-#             num_holes <- length(holes_id)
-#             holes <- matrix(0, nrow=num_holes, ncol=2)
-#             if(num_holes > 0){
-#               for(i in 1:num_holes){
-#                 path  <- t(domain$geometry$edges[domain$geometry$edges_ring==holes_id[i],])[1,]
-#                 path  <- c(path,path[1]) 
-#                 nodes <- domain$geometry$nodes[path,]
-#                 holes[i,] <- st_point_on_surface(st_polygon(list(nodes)))
-#               }
-#             }else{
-#               holes <- NA
-#             }
-#             pslg <- pslg(P = domain$geometry$nodes, PB = domain$geometry$nodes_group,
-#                          S = domain$geometry$edges, SB = domain$geometry$edges_group,
-#                          H = holes) 
-#             triangulation <- triangulate(p = pslg, a = maximum_area, q = minimum_angle)
-#             res <- Mesh(triangulation)
-#             res$geometry <- domain$geometry
-#             res$time_interval <- domain$time_interval
-#             res$crs <- domain$crs
-#             return(res)
-# })
 
 #' create spatio-temporal domain
 #'
diff --git a/R/function_space.R b/R/function_space.R
index 86bedaf..623f913 100644
--- a/R/function_space.R
+++ b/R/function_space.R
@@ -141,7 +141,16 @@ plot.FunctionObject <- function(x, ...){
                             Z=rep(0, times=length(times)),
                             coeff=as.vector(x$coeff[1: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))
+    coeff <- matrix(nrow=nrow(x$FunctionSpace$mesh$get_elements()), 
+                    ncol=length(times))
+    for(t in 1:length(times)){
+      coeff[,t] <- apply(mesh$get_elements(), MARGIN=1, FUN=
+                           function(edge){
+                             mean(x$coeff[edge,t])
+                           })
+    }
+    limits = c(min(coeff), max(coeff))
+    cmin = limits[1]; cmax=limits[2]
     I=x$FunctionSpace$mesh$get_elements()[,1]-1
     J=x$FunctionSpace$mesh$get_elements()[,2]-1
     K=x$FunctionSpace$mesh$get_elements()[,3]-1
diff --git a/R/sf.R b/R/sf.R
index ad55668..2956c83 100644
--- a/R/sf.R
+++ b/R/sf.R
@@ -11,146 +11,116 @@
 #' @importFrom sf st_as_sf
 #' @importFrom sf st_polygon
 #' @importFrom sf st_as_sfc
-#' @name sf
-#' @export
-st_as_sfc.DomainObject <- function(x, ...){
-  polygon_list <- list(mode="list", length=length(x$geometry))
-  for(sub_id in 1:length(x$geometry)){
-    groups_id <- unique(x$geometry[[sub_id]]$edges_ring)
-    n_physical_lines <- length(groups_id)
-    path_list = vector(mode="list", length=n_physical_lines)
-    for(i in 1:n_physical_lines){
-      path  <- t(x$geometry[[sub_id]]$edges[x$geometry[[sub_id]]$edges_ring==groups_id[i],])[1,]
-      path  <- c(path,path[1]) 
-      nodes <- as.matrix(x$coords[path,1:2]) #x$geometry[[sub_id]]$nodes[path,]
-      path_list[[i]] <- nodes
-    }
-    polygon_list[[sub_id]] <- st_cast(st_polygon(path_list),
-                                      to="MULTIPOLYGON")
-  }
-  crs <- x$crs
-  if(is.na(crs)) crs <- NA_crs_
-  
-  if(length(x$geometry)==1){
-    result <- st_sfc(polygon_list[[1]], crs=crs)
-  }else{
-    result <- st_sfc(polygon_list, crs=crs)
-  }
-  return(result)
-}
-
-
-#' @name sf
-#' @importFrom sf st_as_sf
-#' @importFrom sf st_polygon
-#' @importFrom sf st_linestring
-#' @importFrom sf st_point
-#' @importFrom sf st_sfc
+#' @rdname sf
 #' @export
 st_as_sf.DomainObject <- function(x, ...){
   geom_sf <- list()
   for(sub_id in 1:length(x$geometry)){
-    label <- vector(mode="character")
-    sub <- vector(mode="integer")
-    local_id <- vector(mode="integer")
-    group <- vector(mode="numeric")
-    boundary <- vector(mode="integer")
-    bc = vector(mode="character")
-    path_list = vector(mode="list")
-    edge_list = vector(mode="list")
-    pts_list = vector(mode="list")
     
     geometry <- x$geometry[[sub_id]]
     path_id <- unique(geometry$edges_ring)
-    n_physical_lines <- length(path_id)
     
-    path_list_sub = vector(mode="list", length=n_physical_lines)
-    for(i in 1:n_physical_lines){
-      path  <- t(geometry$edges[geometry$edges_ring==path_id[i],][,1]) # )[1,]
+    path_list = vector(mode="list", length=length(path_id))
+    for(i in 1:length(path_id)){
+      path  <- t(geometry$edges[geometry$edges_ring==path_id[i],][,1]) 
       path  <- c(path,path[1]) 
-      nodes <- as.matrix(x$coords[path,1:2]) #geometry$nodes[path,]
-      path_list_sub[[i]] <- st_linestring(nodes)
-    }
-    edge_list_sub = vector(mode="list", length=n_physical_lines)
-    for(i in 1:nrow(geometry$edges)){
-      edge  <- geometry$edges[i,]
-      nodes <- as.matrix(x$coords[edge,1:2]) #geometry$nodes[path,]
-      edge_list_sub[[i]] <- st_linestring(nodes)
+      nodes <- as.matrix(x$coords[path,1:2]) 
+      path_list[[i]] <- st_linestring(nodes)
     }
     
-    groups_pts_id <- unique(geometry$nodes_group)
-    n_physical_pts <- length(groups_pts_id)
-    pts_list_sub = vector(mode="list", length=n_physical_pts)
-    for(i in 1:nrow(geometry$nodes)){
-      node  <- geometry$nodes[i,]
-      pts_list_sub[[i]] <- st_point(t(as.matrix(node)))
-    }
+    # edges to st_linestring ! (FAST)
+    edge_list <- lapply(as.list(as.data.frame(t(geometry$edges))), FUN=
+                          function(edge){
+                            st_linestring(as.matrix(x$coords[edge,1:2]))
+                          })
     
-    label= append(label, c(rep("path", times=n_physical_lines),
-                           rep("edge", times=nrow(geometry$edges)),
-                           rep("node", times=nrow(geometry$nodes))))
-    group = append(group, c(path_id , geometry$edges_group, geometry$nodes_group))
-    local_id = append(local_id, 
-                      c(1:n_physical_lines,1:nrow(geometry$edges), 1:nrow(geometry$nodes)))
+    # nodes to st_point ! (FAST)
+    pts_list <- lapply(as.list(as.data.frame(t(geometry$nodes))), st_point)
     
-    path_list = append(path_list, path_list_sub) 
-    edge_list = append(edge_list, edge_list_sub)
-    pts_list  = append(pts_list, pts_list_sub)
+    label= c(rep("path", times=length(path_id)),
+             rep("edge", times=length(edge_list)),
+             rep("node", times=length(pts_list)))
+    group = c(path_id , geometry$edges_group, geometry$nodes_group)
+    id =  c(1:length(path_id),
+            1:nrow(geometry$edges), 
+            1:nrow(geometry$nodes))
     
     group <- as.factor(group)
     crs <- x$crs
     if(is.na(crs)) crs <- NA_crs_
     geom_sfc <- st_sfc( append(append(path_list, edge_list), pts_list), crs=crs)
     geom_sf[[sub_id]] <- st_as_sf(data.frame(label= label,
-                                             group = group, local_id=local_id),
+                                             group = group, id=id),
                                   geometry = geom_sfc)
     
   }
   if(length(x$geometry)==1) geom_sf <- geom_sf[[1]]
   return(geom_sf)
 }
-# st_as_sf.DomainObject <- function(x, ...){
-#   groups_id <- unique(x$geometry$edges_ring)
-#   n_physical_lines <- length(groups_id)
-#   path_list = vector(mode="list", length=n_physical_lines)
-#   for(i in 1:n_physical_lines){
-#     path  <- t(x$geometry$edges[x$geometry$edges_ring==groups_id[i],])[1,]
-#     path  <- c(path,path[1]) 
-#     nodes <- x$geometry$nodes[path,]
-#     path_list[[i]] <- nodes
-#   }
-#   st_sfc( st_polygon(path_list) )
-# }
 
 #' @importFrom sf st_as_sf
 #' @importFrom sf st_polygon
+#' @importFrom sf st_linestring
+#' @importFrom sf st_point
 #' @importFrom sf st_sfc
-#' @name sf
+#' @rdname sf
 #' @export
-st_as_sfc.MeshObject <- function(x, ...){
-  polygon_list <- list()
-  for(e in 1:nrow(x$get_elements())){
-    element_sf <- st_polygon(list(rbind(x$get_nodes()[x$get_elements()[e,],],
-                                        x$get_nodes()[x$get_elements()[e,1],])))
-    polygon_list[[e]]  <- element_sf
+st_as_sfc.DomainObject <- function(x, ...){
+  polygon_list <- list(mode="list", length=length(x$geometry))
+  for(sub_id in 1:length(x$geometry)){
+    geometry <- x$geometry[[sub_id]]
     
+    path_id <- unique(geometry$edges_ring)
+    path_list = vector(mode="list", length=length(path_id))
+    for(i in 1:length(path_id)){
+      path  <- t(geometry$edges[geometry$edges_ring==path_id[i],])[1,]
+      path  <- c(path,path[1]) 
+      nodes <- as.matrix(x$coords[path,1:2]) 
+      path_list[[i]] <- nodes
+    }
+    polygon_list[[sub_id]] <- st_cast(st_polygon(path_list),
+                                      to="MULTIPOLYGON")
+  }
+  crs <- x$crs
+  if(is.na(crs)) crs <- NA_crs_
+  
+  if(length(x$geometry)==1){
+    result <- st_sfc(polygon_list[[1]], crs=crs)
+  }else{
+    result <- st_sfc(polygon_list, crs=crs)
   }
+  return(result)
+}
+
+#' @importFrom sf st_as_sf
+#' @importFrom sf st_polygon
+#' @importFrom sf st_sfc
+#' @rdname sf
+#' @export
+st_as_sfc.MeshObject <- function(x, ...){
+
+  polygon_list <- apply(x$get_elements(), MARGIN=1, FUN=function(x){
+    st_cast(st_linestring(mesh$get_nodes()[x,]),
+            to ="POLYGON"
+    )
+  })
+  
   crs <- x$crs
   if(is.na(crs)) crs <- NA_crs_
   mesh_sf <- st_sfc(polygon_list, crs = crs)
   mesh_sf
 }
 
-#' @name sf
 #' @importFrom sf st_crs
+#' @rdname sf
 #' @export
 st_crs.DomainObject <- function(x, ...){
   #st_crs(x$crs, ...)
   st_crs(st_as_sfc(x), ...)
 }
 
-#' @name sf
 #' @importFrom sf st_crs<- st_crs
+#' @rdname sf
 #' @export
 `st_crs<-.DomainObject` = function(x, value) {
   #st_crs(x$crs) = value
diff --git a/tests/build_domain.R b/tests/build_domain.R
index df3f1c5..4a21c9b 100644
--- a/tests/build_domain.R
+++ b/tests/build_domain.R
@@ -29,12 +29,12 @@ st_crs(domain_sf)
 domain_sf <- st_as_sf(domain) # + dataset
 plot(domain_sf)
 domain_sf %>% filter(label=="edge") %>% 
-              select(local_id)
+              select(id)
 plot(domain_sf %>% filter(label=="edge") %>% 
-       select(local_id), lwd=3)
+       select(id), lwd=3)
 
-plot(domain_sf %>% filter(label=="edge" & local_id >3) %>% 
-       select(local_id), lwd=3)
+plot(domain_sf %>% filter(label=="edge" & id >3) %>% 
+       select(id), lwd=3)
 #plot(domain_sf %>% select(bc) %>% filter(!is.na(bc)))
 #plot(domain_sf %>% select() %>% filter(!is.na(bc)))
 
@@ -44,7 +44,7 @@ mesh_sf <- st_as_sfc(mesh)
 plot(mesh_sf, col="red")
 
 # espolare st_triangulate / st_triangulate_constrained
-plot( st_triangulate(mesh_sf), col="red" )
+#plot( st_triangulate(mesh_sf), col="red" )
 
 # 2.
 nodes_ext <- rbind(c(0, 0), c(1, 0), c(1, 1), c(0, 1))