diff --git a/Manifest.toml b/Manifest.toml
index f3074c911..135c36e16 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -1,5 +1,11 @@
 # This file is machine-generated - editing it directly is not advised
 
+[[ArrayInterface]]
+deps = ["LinearAlgebra", "Requires", "SparseArrays"]
+git-tree-sha1 = "981354dab938901c2b607a213e62d9defa50b698"
+uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+version = "1.2.1"
+
 [[Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
@@ -15,6 +21,12 @@ git-tree-sha1 = "c7361ce8a2129f20b0e05a89f7070820cfed6648"
 uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
 version = "0.5.6"
 
+[[Calculus]]
+deps = ["Compat"]
+git-tree-sha1 = "bd8bbd105ba583a42385bd6dc4a20dad8ab3dc11"
+uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
+version = "0.5.0"
+
 [[CodecZlib]]
 deps = ["BinaryProvider", "Libdl", "TranscodingStreams"]
 git-tree-sha1 = "05916673a2627dd91b4969ff8ba6941bc85a960e"
@@ -27,6 +39,12 @@ git-tree-sha1 = "50b3ae4d643dc27eaff69fb6be06ee094d5500c9"
 uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 version = "0.7.0"
 
+[[CommonSubexpressions]]
+deps = ["Test"]
+git-tree-sha1 = "efdaf19ab11c7889334ca247ff4c9f7c322817b0"
+uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
+version = "0.2.0"
+
 [[Compat]]
 deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
 git-tree-sha1 = "84aa74986c5b9b898b0d1acaf3258741ee64754f"
@@ -47,6 +65,30 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 deps = ["Mmap"]
 uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 
+[[DiffEqDiffTools]]
+deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"]
+git-tree-sha1 = "21b855cb29ec4594f9651e0e9bdc0cdcfdcd52c1"
+uuid = "01453d9d-ee7c-5054-8395-0335cb756afa"
+version = "1.3.0"
+
+[[DiffResults]]
+deps = ["Compat", "StaticArrays"]
+git-tree-sha1 = "34a4a1e8be7bc99bc9c611b895b5baf37a80584c"
+uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
+version = "0.0.4"
+
+[[DiffRules]]
+deps = ["Random", "Test"]
+git-tree-sha1 = "dc0869fb2f5b23466b32ea799bd82c76480167f7"
+uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
+version = "0.0.10"
+
+[[Distances]]
+deps = ["LinearAlgebra", "Statistics"]
+git-tree-sha1 = "23717536c81b63e250f682b0e0933769eecd1411"
+uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+version = "0.8.2"
+
 [[Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
@@ -57,6 +99,12 @@ git-tree-sha1 = "769ac4057ed875f433372e9a571a2cb86347d1be"
 uuid = "442a2c76-b920-505d-bb47-c5924d526838"
 version = "0.3.3"
 
+[[ForwardDiff]]
+deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "InteractiveUtils", "LinearAlgebra", "NaNMath", "Random", "SparseArrays", "SpecialFunctions", "StaticArrays", "Test"]
+git-tree-sha1 = "4c4d727f1b7e0092134fabfab6396b8945c1ea5b"
+uuid = "f6369f11-7733-5829-9624-2563aa707210"
+version = "0.10.3"
+
 [[InteractiveUtils]]
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
@@ -79,6 +127,12 @@ git-tree-sha1 = "aeec7a341652d47bc773475a42952fa78eccd7cc"
 uuid = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
 version = "0.8.0"
 
+[[LineSearches]]
+deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf", "Test"]
+git-tree-sha1 = "54eb90e8dbe745d617c78dee1d6ae95c7f6f5779"
+uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
+version = "7.0.1"
+
 [[LinearAlgebra]]
 deps = ["Libdl"]
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -93,12 +147,36 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 [[Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
+[[NLSolversBase]]
+deps = ["Calculus", "DiffEqDiffTools", "DiffResults", "Distributed", "ForwardDiff"]
+git-tree-sha1 = "c430bd3f2dfcffc30688cf4a9cb61535e8d85f65"
+uuid = "d41bc354-129a-5804-8e4c-c37616107c6c"
+version = "7.4.1"
+
+[[NLsolve]]
+deps = ["DiffEqDiffTools", "Distances", "ForwardDiff", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"]
+git-tree-sha1 = "c9578878f56f425a2160f5b436c7f79a154d154c"
+uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
+version = "4.1.0"
+
+[[NaNMath]]
+deps = ["Compat"]
+git-tree-sha1 = "ce3b85e484a5d4c71dd5316215069311135fa9f2"
+uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+version = "0.3.2"
+
 [[OrderedCollections]]
 deps = ["Random", "Serialization", "Test"]
 git-tree-sha1 = "c4c13474d23c60d20a67b217f1d7f22a40edf8f1"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 version = "1.1.0"
 
+[[Parameters]]
+deps = ["OrderedCollections"]
+git-tree-sha1 = "1dfd7cd50a8eb06ef693a4c2bbe945943cd000c5"
+uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+version = "0.11.0"
+
 [[Parsers]]
 deps = ["Dates", "Test"]
 git-tree-sha1 = "ef0af6c8601db18c282d092ccbd2f01f3f0cd70b"
@@ -139,6 +217,12 @@ git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"
 uuid = "189a3867-3050-52da-a836-e630ba90ab69"
 version = "0.2.0"
 
+[[Requires]]
+deps = ["Test"]
+git-tree-sha1 = "f6fbf4ba64d295e146e49e021207993b6b48c7d1"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "0.5.2"
+
 [[SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 
diff --git a/NEWS.md b/NEWS.md
index 214a704c8..aee698b76 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -7,6 +7,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
+- Now, the `tags` argument is optional when constucting `SkeletonTriangulation` and `BoundaryTriangulation` objects from a `DiscreteModel`. Since commit [e6424a3](https://github.com/gridap/Gridap.jl/commit/e6424a304feb38547241e86de07a821e26344a7e).
+- Added `mean` operator for quantities restricted to a `SkeletonTriangulation`. Since commit [83798b4](https://github.com/gridap/Gridap.jl/commit/83798b4f38aaf482b968ffd0359eb75c79a21385).  
+- Extended `NormalVector` to `SkeletonTriangulations`. Since commit [5fb8487](https://github.com/gridap/Gridap.jl/commit/5fb84871128c4388559cc5052d9ff00f0be19462). 
+- Now, `TrialFESpaces` can be constructed from values instead of functions if the corresponding Dirichlet conditions are constant. Since commit [bae237e](https://github.com/gridap/Gridap.jl/commit/bae237e881db6569622f3559f82bcc3999560526).
+- Added the possibility of adding new tags to a `FaceLabels` object via the function `add_tag_from_tags!` and using it to construct FE spaces. Since commit [e9dfac4](https://github.com/gridap/Gridap.jl/commit/e9dfac4489047c0b7e1c62507f4335e9fc76dfd8).
+- Added `BackslashSolver` to facilitate the usage in Gridap of the build-in Julia backslash linear solver. Since commit [8e3a9b7](https://github.com/gridap/Gridap.jl/commit/8e3a9b71c64b032c5a572a7ef696f4cbf875190b).
+- Added `JuliaNLSolver` to facilitat the usage in Gridap of the non-linear solvers available in the official Julia package `NLsolve`. Since commit [e5a933f](https://github.com/gridap/Gridap.jl/commit/e5a933f3093faea221a50bdd796d7f02113ed52c).
+
 ### Changed
 ### Removed
 ### Deprecated
diff --git a/Project.toml b/Project.toml
index 5fa52890d..49e37ee62 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
diff --git a/src/Algebra/JuliaNLSolvers.jl b/src/Algebra/JuliaNLSolvers.jl
new file mode 100644
index 000000000..32a8bb567
--- /dev/null
+++ b/src/Algebra/JuliaNLSolvers.jl
@@ -0,0 +1,66 @@
+module JuliaNLSolvers
+
+using Gridap
+using NLsolve
+
+export JuliaNLSolver
+
+import Gridap: solve!
+
+mutable struct JuliaNLSolver <: NonLinearSolver
+  ls::LinearSolver
+  kwargs::Dict
+end
+
+function JuliaNLSolver(ls::LinearSolver;kwargs...)
+  @assert ! haskey(kwargs,:linsolve) "linsolve cannot be used here. It is managed internally"
+  JuliaNLSolver(ls,kwargs)
+end
+
+function JuliaNLSolver(;kwargs...)
+  ls = BackslashSolver()
+  JuliaNLSolver(ls;kwargs...)
+end
+
+mutable struct JuliaNLSolversCache
+  df::OnceDifferentiable
+  ss::SymbolicSetup
+  result
+end
+
+function solve!(x::AbstractVector,nls::JuliaNLSolver,op::NonLinearOperator)
+  cache = _setup_cache(x,nls,op)
+  solve!(x,nls,op,cache)
+  cache
+end
+
+function solve!(
+  x::AbstractVector,nls::JuliaNLSolver,op::NonLinearOperator,cache::JuliaNLSolversCache)
+  df = cache.df
+  ss = cache.ss
+  kwargs = nls.kwargs
+  function linsolve!(x,A,b)
+    ns = numerical_setup(ss,A)
+    solve!(x,ns,A,b)
+  end
+  r = nlsolve(df,x;linsolve=linsolve!,kwargs...)
+  cache.result = r
+  x[:] .= r.zero
+end
+
+function _setup_cache(x0,nls,op)
+
+  f!(r,x) = residual!(r,op,x)
+  j!(j,x) = jacobian!(j,op,x)
+  
+  f0 = residual(op,x0)
+  j0 = jacobian(op,x0)
+  
+  df = OnceDifferentiable(f!,j!,x0,f0,j0)
+
+  ss = symbolic_setup(nls.ls,j0)
+
+  JuliaNLSolversCache(df,ss,nothing)
+end
+
+end # module
diff --git a/src/Algebra/LinearSolvers.jl b/src/Algebra/LinearSolvers.jl
index 10db58d2c..159ad4911 100644
--- a/src/Algebra/LinearSolvers.jl
+++ b/src/Algebra/LinearSolvers.jl
@@ -15,6 +15,7 @@ export numerical_setup
 export numerical_setup!
 export test_linear_solver
 export LUSolver
+export BackslashSolver
 
 abstract type LinearSolver end
 
@@ -80,4 +81,26 @@ function solve!(
   x .= y
 end
 
+"""
+Wrapper of the backslash solver available in julia
+This is typically faster than LU for a single solve
+"""
+struct BackslashSolver <: LinearSolver end
+
+struct BackslashSymbolicSetup <: SymbolicSetup end
+
+struct BackslashNumericalSetup <: NumericalSetup end
+
+symbolic_setup(::BackslashSolver,mat::AbstractMatrix) = BackslashSymbolicSetup()
+
+numerical_setup(::BackslashSymbolicSetup,mat::AbstractMatrix) = BackslashNumericalSetup()
+
+function numerical_setup!(ns::BackslashNumericalSetup, mat::AbstractMatrix)
+end
+
+function solve!(
+  x::AbstractVector,ns::BackslashNumericalSetup,A::AbstractMatrix,b::AbstractVector)
+   copyto!(x, A\b)
+end
+
 end
diff --git a/src/Algebra/files.jl b/src/Algebra/files.jl
index ae764a96f..97fea0564 100644
--- a/src/Algebra/files.jl
+++ b/src/Algebra/files.jl
@@ -4,3 +4,6 @@ include("LinearSolvers.jl")
 
 include("NonLinearSolvers.jl")
 @reexport using Gridap.NonLinearSolvers
+
+include("JuliaNLSolvers.jl")
+@reexport using Gridap.JuliaNLSolvers
diff --git a/src/FESpaces/CLagrangianFESpaces.jl b/src/FESpaces/CLagrangianFESpaces.jl
index 1cf2d7a00..5238671af 100644
--- a/src/FESpaces/CLagrangianFESpaces.jl
+++ b/src/FESpaces/CLagrangianFESpaces.jl
@@ -81,12 +81,26 @@ function CLagrangianFESpace(
 end
 
 function CLagrangianFESpace(
-  ::Type{T},model::DiscreteModel,order,diritags=Int[],dirimasks=nothing) where T
+  ::Type{T},
+  model::DiscreteModel,
+  order::Integer,
+  diritags=Int[],
+  dirimasks=nothing) where T
+  CLagrangianFESpace(T,model,FaceLabels(model),order,diritags,dirimasks)
+end
+
+function CLagrangianFESpace(
+  ::Type{T},
+  model::DiscreteModel,
+  facelabels::FaceLabels,
+  order::Integer,
+  diritags=Int[],
+  dirimasks=nothing) where T
 
-  _diri_tags = _setup_tags(model,diritags)
+  _diri_tags = _setup_tags(facelabels,diritags)
   _diri_masks = _setup_masks(T,_diri_tags,dirimasks)
 
-  grid, node_to_label, tag_to_labels = _setup_grid(model,order)
+  grid, node_to_label, tag_to_labels = _setup_grid(model,facelabels,order)
 
   CLagrangianFESpace(
     T,grid,node_to_label,tag_to_labels,_diri_tags,_diri_masks)
@@ -111,6 +125,14 @@ apply_constraints_cols(
 celldofids(fesp::CLagrangianFESpace) = fesp.cell_to_dofs
 
 function interpolate_values(fesp::CLagrangianFESpace,fun::Function)
+  _interpolate_values(fesp,fun)
+end
+
+function interpolate_values(fesp::CLagrangianFESpace{D,Z,T},val::T) where {D,Z,T}
+  _interpolate_values(fesp,val)
+end
+
+function _interpolate_values(fesp,fun)
   zh = zero(fesp)
   fdof_to_val = free_dofs(zh)
   ddof_to_val = diri_dofs(zh)
@@ -119,6 +141,14 @@ function interpolate_values(fesp::CLagrangianFESpace,fun::Function)
 end
 
 function interpolate_diri_values(fesp::CLagrangianFESpace, funs::Vector{<:Function})
+  _interpolate_diri_values(fesp, funs)
+end
+
+function interpolate_diri_values(fesp::CLagrangianFESpace{D,Z,T},vals::Vector{T}) where {D,Z,T}
+  _interpolate_diri_values(fesp, vals)
+end
+
+function _interpolate_diri_values(fesp, funs)
 
   T = value_type(fesp)
   E = eltype(T)
@@ -359,7 +389,7 @@ end
 function _fill_interpolated_vals!(fdof_to_val,ddof_to_val,fesp,fun)
 
   x = points(fesp.grid)
-  fx = fun.(x)
+  fx = _eval_bc(fun,x)
   T = value_type(fesp)
   @assert eltype(fx) == T
   z = zero(T)
@@ -374,6 +404,10 @@ function _fill_interpolated_vals!(fdof_to_val,ddof_to_val,fesp,fun)
 
 end
 
+_eval_bc(fun::Function,x) = fun.(x)
+
+_eval_bc(val,x) = fill(val,length(x))
+
 function  _fill_interpolated_vals_kernel!(
   fdof_to_val,
   ddof_to_val,
@@ -418,7 +452,7 @@ function _fill_diri_values!(
     label = node_to_label[node]
     if label in tag_to_labels[diritag]
       x = node_to_coords[node]
-      fx = fun(x)
+      fx = _eval(fun,x)
       comp_to_dof = node_and_comp_to_dof[node]
       for comp in eachindex(fx)
         mask = dirimask[comp]
@@ -435,6 +469,10 @@ function _fill_diri_values!(
 
 end
 
+@inline _eval(fun::Function,x) = fun(x)
+
+@inline _eval(val,x) = val
+
 function _same_order(order,co)
   @notimplemented
 end
@@ -443,9 +481,8 @@ function _same_order(order,co::ConstantCellValue)
   order == co.value
 end
 
-function _setup_grid(model,order)
+function _setup_grid(model,facelabels,order)
 
-  facelabels = FaceLabels(model)
   _grid = Grid(model)
 
   if _same_order(order,cellorders(_grid))
diff --git a/src/FESpaces/ConformingFESpaces.jl b/src/FESpaces/ConformingFESpaces.jl
index 392662b05..ad4737108 100644
--- a/src/FESpaces/ConformingFESpaces.jl
+++ b/src/FESpaces/ConformingFESpaces.jl
@@ -58,15 +58,27 @@ function ConformingFESpace(
   return ConformingFESpace(reffe, trian, graph, labels, ())
 end
 
-function ConformingFESpace(::Type{T},model::DiscreteModel{D},order,diri_tags) where {D,T}
+function ConformingFESpace(
+  ::Type{T}, model::DiscreteModel{D}, order::Integer, diri_tags) where {D,T}
+
+  labels = FaceLabels(model)
+  ConformingFESpace(T,model,labels,order,diri_tags)
+end
+
+function ConformingFESpace(
+  ::Type{T},
+  model::DiscreteModel{D},
+  labels::FaceLabels,
+  order::Integer,
+  diri_tags) where {D,T}
+
   grid = Grid(model,D)
   trian = Triangulation(grid)
   graph = GridGraph(model)
-  labels = FaceLabels(model)
   orders = fill(order,D)
   polytope = _polytope(celltypes(grid))
   fe = LagrangianRefFE(T,polytope, orders)
-  _diri_tags = _setup_tags(model,diri_tags)
+  _diri_tags = _setup_tags(labels,diri_tags)
   ConformingFESpace(fe,trian,graph,labels,_diri_tags)
 end
 
@@ -99,10 +111,19 @@ function interpolate_values(this::ConformingFESpace,f::Function)
   _interpolate_values(this,f)
 end
 
+function interpolate_values(this::ConformingFESpace{D,Z,T},val::T) where {D,Z,T}
+  fun(x) = val
+  interpolate_values(this,fun)
+end
+
 function interpolate_diri_values(this::ConformingFESpace, funs::Vector{<:Function})
   _interpolate_diri_values(this,funs)
 end
 
+function interpolate_diri_values(this::ConformingFESpace{D,Z,T}, vals::Vector{T}) where {D,Z,T}
+  _interpolate_diri_values(this,vals)
+end
+
 function CellField(
   fespace::ConformingFESpace{D,Z,T},
   free_dofs::AbstractVector{E},
diff --git a/src/FESpaces/ConstrainedFESpaces.jl b/src/FESpaces/ConstrainedFESpaces.jl
index 60e9bf964..1021da8d4 100644
--- a/src/FESpaces/ConstrainedFESpaces.jl
+++ b/src/FESpaces/ConstrainedFESpaces.jl
@@ -85,10 +85,19 @@ function interpolate_values(this::ConstrainedFESpace,fun::Function)
   (freevals[this.is_free], dirivals)
 end
 
+function interpolate_values(this::ConstrainedFESpace{D,Z,T},val::T) where {D,Z,T}
+  freevals, dirivals = interpolate_values(this.fespace,val)
+  (freevals[this.is_free], dirivals)
+end
+
 function interpolate_diri_values(this::ConstrainedFESpace, funs::Vector{<:Function})
   interpolate_diri_values(this.fespace,funs)
 end
 
+function interpolate_diri_values(this::ConstrainedFESpace{D,Z,T}, vals::Vector{T}) where {D,Z,T}
+  interpolate_diri_values(this.fespace,vals)
+end
+
 function CellField(
   this::ConstrainedFESpace{D,Z,T},
   free_dofs::AbstractVector{E},
diff --git a/src/FESpaces/DLagrangianFESpaces.jl b/src/FESpaces/DLagrangianFESpaces.jl
index 64fa35dbd..588cda39c 100644
--- a/src/FESpaces/DLagrangianFESpaces.jl
+++ b/src/FESpaces/DLagrangianFESpaces.jl
@@ -9,6 +9,8 @@ using Gridap.CLagrangianFESpaces: _S
 using Gridap.CLagrangianFESpaces: _compute_comp_to_dof
 using Gridap.CLagrangianFESpaces: _setup_grid
 using Gridap.CLagrangianFESpaces: _setup_masks
+using Gridap.CLagrangianFESpaces: _eval
+using Gridap.CLagrangianFESpaces: _eval_bc
 using Gridap.ConformingFESpaces: _CellField
 using Gridap.BoundaryGrids: _setup_tags
 
@@ -80,12 +82,26 @@ function DLagrangianFESpace(
 end
 
 function DLagrangianFESpace(
-  ::Type{T},model::DiscreteModel,order,diritags=Int[],dirimasks=nothing) where T
+  ::Type{T},
+  model::DiscreteModel,
+  order::Integer,
+  diritags=Int[],
+  dirimasks=nothing) where T
+  DLagrangianFESpace(T,model,FaceLabels(model),order,diritags,dirimasks)
+end
 
-  _diri_tags = _setup_tags(model,diritags)
+function DLagrangianFESpace(
+  ::Type{T},
+  model::DiscreteModel,
+  facelabels::FaceLabels,
+  order::Integer,
+  diritags=Int[],
+  dirimasks=nothing) where T
+
+  _diri_tags = _setup_tags(facelabels,diritags)
   _diri_masks = _setup_masks(T,_diri_tags,dirimasks)
 
-  grid, node_to_label, tag_to_labels = _setup_grid(model,order)
+  grid, node_to_label, tag_to_labels = _setup_grid(model,facelabels,order)
 
   DLagrangianFESpace(
     T,grid,node_to_label,tag_to_labels,_diri_tags,_diri_masks)
@@ -118,6 +134,14 @@ function CellField(
 end
 
 function interpolate_values(fesp::DLagrangianFESpace,fun::Function)
+  _interpolate_values(fesp,fun)
+end
+
+function interpolate_values(fesp::DLagrangianFESpace{D,Z,T},val::T) where {D,Z,T}
+  _interpolate_values(fesp,val)
+end
+
+function _interpolate_values(fesp::DLagrangianFESpace,fun)
 
   zh = zero(fesp)
   fdof_to_val = free_dofs(zh)
@@ -126,7 +150,7 @@ function interpolate_values(fesp::DLagrangianFESpace,fun::Function)
   lnode_and_comp_to_ldof = fesp.reffe.dofbasis.node_and_comp_to_dof
   cell_to_dofs = fesp.cell_to_dofs
   cell_to_nodes = cells(fesp.grid)
-  node_comp_to_val = fun.(node_to_coords)
+  node_comp_to_val = _eval_bc(fun,node_to_coords)
 
   _fill_interpolated_vals!(
     fdof_to_val,
@@ -140,6 +164,14 @@ function interpolate_values(fesp::DLagrangianFESpace,fun::Function)
 end
 
 function interpolate_diri_values(fesp::DLagrangianFESpace, funs::Vector{<:Function})
+  _interpolate_diri_values(fesp, funs)
+end
+
+function interpolate_diri_values(fesp::DLagrangianFESpace{D,Z,T}, vals::Vector{T}) where {D,Z,T}
+  _interpolate_diri_values(fesp, vals)
+end
+
+function _interpolate_diri_values(fesp::DLagrangianFESpace, funs)
 
   T = value_type(fesp)
   E = eltype(T)
@@ -202,7 +234,7 @@ function _fill_diri_values!(
       if label in tag_to_labels[diritag]
 
         coords = node_to_coords[node]
-        comp_to_val = fun(coords)
+        comp_to_val = _eval(fun,coords)
 
         for comp in eachindex(dirimask)
           mask = dirimask[comp]
diff --git a/src/FESpaces/DiscFESpaces.jl b/src/FESpaces/DiscFESpaces.jl
index 275c3c81c..125b91552 100644
--- a/src/FESpaces/DiscFESpaces.jl
+++ b/src/FESpaces/DiscFESpaces.jl
@@ -71,7 +71,20 @@ function interpolate_values(this::DiscFESpace,fun::Function)
   _setup_interpolate_values(this, fun)
 end
 
+function interpolate_values(this::DiscFESpace{D,Z,T},val::T) where {D,Z,T}
+  fun(x) = val
+  _setup_interpolate_values(this, fun)
+end
+
 function interpolate_diri_values(this::DiscFESpace, funs::Vector{<:Function})
+  _interpolate_diri_values(this::DiscFESpace, funs)
+end
+
+function interpolate_diri_values(this::DiscFESpace{D,Z,T}, vals::Vector{T}) where {D,Z,T}
+  _interpolate_diri_values(this::DiscFESpace, vals)
+end
+
+function _interpolate_diri_values(this::DiscFESpace, funs)
   @assert length(funs) == 0
   T = value_type(this)
   E = dof_type(T)
diff --git a/src/FESpaces/FEBases.jl b/src/FESpaces/FEBases.jl
index 4a7201207..b29bf10e7 100644
--- a/src/FESpaces/FEBases.jl
+++ b/src/FESpaces/FEBases.jl
@@ -11,6 +11,7 @@ import Base: div
 import Gridap: trace
 import Gridap: curl
 import Gridap: inner
+import Gridap: varinner
 import Base: +, -, *
 import Gridap: restrict
 import Gridap: Triangulation
@@ -60,6 +61,10 @@ function inner(a::FEBasis,b::CellField)
   varinner(a.cellbasis,b)
 end
 
+function varinner(a::FEBasis,b::CellField)
+  varinner(a.cellbasis,b)
+end
+
 function inner(a::FEBasis,f::Function)
   b = CellField(a.trian,f)
   inner(a,b)
@@ -69,6 +74,10 @@ function inner(a::FEBasis,b::FEBasis)
   varinner(a.cellbasis,b.cellbasis)
 end
 
+function varinner(a::FEBasis,b::FEBasis)
+  varinner(a.cellbasis,b.cellbasis)
+end
+
 Triangulation(a::FEBasis) = a.trian
 
 function CellBasis(
@@ -86,7 +95,14 @@ function restrict(feb::FEBasis,trian::BoundaryTriangulation)
 end
 
 function restrict(feb::FEBasis,trian::SkeletonTriangulation)
-  restrict(feb.cellbasis,trian)
+  sp = restrict(feb.cellbasis,trian)
+  _new_sp(sp,trian)
+end
+
+function _new_sp(sp::SkeletonPair{Z,T,N},trian) where {Z,T,N}
+  b1 = FEBasis(sp.cellfield1,trian)
+  b2 = FEBasis(sp.cellfield2,trian)
+  SkeletonPair{Z,T,N}(b1,b2)
 end
 
 end # module
diff --git a/src/FESpaces/FESpaces.jl b/src/FESpaces/FESpaces.jl
index 39f8bad73..2bbf934b4 100644
--- a/src/FESpaces/FESpaces.jl
+++ b/src/FESpaces/FESpaces.jl
@@ -71,6 +71,10 @@ function interpolate_values(::FESpace,::Function)::Tuple{Vector{E},Vector{E}} wh
   @abstractmethod
 end
 
+function interpolate_values(::FESpace{D,Z,T},::T) where {D,Z,T}
+  @abstractmethod
+end
+
 """
 Given a FESpace and its array of labels that represent the Dirichlet boundary
 in the geometrical model and an array of functions for every Dirichlet
@@ -81,11 +85,20 @@ function interpolate_diri_values(::FESpace, funs::Vector{<:Function})::Vector{E}
   @abstractmethod
 end
 
+function interpolate_diri_values( ::FESpace{D,Z,T}, ::Vector{T}) where {D,Z,T}
+  @abstractmethod
+end
+
 function interpolate_diri_values(this::FESpace, fun::Function)
   tags = diri_tags(this)
   interpolate_diri_values(this,fill(fun,length(tags)))
 end
 
+function interpolate_diri_values(this::FESpace{D,Z,T}, val::T) where {D,Z,T}
+  tags = diri_tags(this)
+  interpolate_diri_values(this,fill(val,length(tags)))
+end
+
 """
 Returns the CellField that represents the FE function thorough its free and
 dirichlet values. E = eltype(T)
@@ -111,27 +124,40 @@ function TestFESpace(this::FESpace{D,Z,T}) where {D,Z,T}
   return FESpaceWithDirichletData(this, dv)
 end
 
-function TrialFESpace( this::FESpace, funs::Vector{<:Function}=Function[]) where {D}
+function TrialFESpace( this::FESpace, funs::Vector{<:Function}=Function[])
+  _TrialFESpace(this,funs)
+end
+
+function TrialFESpace( this::FESpace{D,Z,T}, vals::Vector{T}) where {D,Z,T}
+  _TrialFESpace(this,vals)
+end
+
+function _TrialFESpace( this, funs)
   tags = diri_tags(this)
   @assert length(tags) == length(funs)
   dv = interpolate_diri_values(this,funs)
   return FESpaceWithDirichletData(this, dv)
 end
 
-function TrialFESpace( this::FESpace, fun::Function) where {D}
+function TrialFESpace( this::FESpace, fun::Function)
   dv = interpolate_diri_values(this,fun)
   return FESpaceWithDirichletData(this, dv)
 end
 
+function TrialFESpace( this::FESpace{D,Z,T}, val::T) where {D,Z,T}
+  dv = interpolate_diri_values(this,val)
+  return FESpaceWithDirichletData(this, dv)
+end
+
 # Testers
 
 function test_fe_space(
-  fespace::FESpace{D,Z},
+  fespace::FESpace{D,Z,T},
   nfree::Integer,
   ndiri::Integer,
   cellmat::CellMatrix,
   cellvec::CellVector,
-  ufun::Function) where {D,Z}
+  ufun::Function) where {D,Z,T}
 
   @test num_free_dofs(fespace) == nfree
   @test num_diri_dofs(fespace) == ndiri
@@ -160,6 +186,13 @@ function test_fe_space(
   @test length(freevals) == nfree
   @test length(dirivals) == ndiri
 
+  z = zero(T)
+  freevals, dirivals = interpolate_values(fespace,z)
+  @test isa(freevals,AbstractVector)
+  @test isa(dirivals,AbstractVector)
+  @test length(freevals) == nfree
+  @test length(dirivals) == ndiri
+
   dirivals = interpolate_diri_values(fespace,ufun)
   @test isa(dirivals,AbstractVector)
   @test length(dirivals) == ndiri
@@ -170,6 +203,10 @@ function test_fe_space(
   trian = Triangulation(fespace)
   @test isa(trian,Triangulation{Z,D})
 
+  dirivals = interpolate_diri_values(fespace,z)
+  @test isa(dirivals,AbstractVector)
+  @test length(dirivals) == ndiri
+
 end
 
 # Helpers
@@ -216,6 +253,11 @@ function interpolate_values(f::FESpaceWithDirichletData,fun::Function)
   free_vals, f.diri_dofs
 end
 
+function interpolate_values(f::FESpaceWithDirichletData{D,Z,T},val::T) where {D,Z,T}
+  free_vals, _ = interpolate_values(f.fespace,val)
+  free_vals, f.diri_dofs
+end
+
 function CellField(
   f::FESpaceWithDirichletData{D,Z,T},free_dofs::AbstractVector{E},diri_dofs::AbstractVector{E})where {D,Z,T,E}
   CellField(f.fespace,free_dofs,f.diri_dofs)
@@ -231,4 +273,8 @@ function interpolate_diri_values(f::FESpaceWithDirichletData, funs::Vector{<:Fun
   f.diri_dofs
 end
 
+function interpolate_diri_values(f::FESpaceWithDirichletData{D,Z,T}, vals::Vector{T}) where {D,Z,T}
+  f.diri_dofs
+end
+
 end # module FESpaces
diff --git a/src/Fields/FieldValues.jl b/src/Fields/FieldValues.jl
index e752a7905..25e0b44c5 100644
--- a/src/Fields/FieldValues.jl
+++ b/src/Fields/FieldValues.jl
@@ -17,6 +17,9 @@ export symmetic_part
 
 import TensorValues: meas
 import Base: adjoint
+import LinearAlgebra: tr
+import LinearAlgebra: dot
+import LinearAlgebra: norm
 
 
 """
@@ -57,7 +60,7 @@ end
   Meta.parse(str[1:(end-1)])
 end
 
-const tr = trace
+@inline tr(v::TensorValue) = trace(v)
 
 @generated function symmetic_part(v::TensorValue{D}) where D
   str = "("
@@ -75,4 +78,8 @@ function adjoint(v::TensorValue)
   TensorValue(t)
 end
 
+@inline dot(u::VectorValue,v::VectorValue) = inner(u,v)
+
+@inline norm(u::VectorValue) = sqrt(inner(u,u))
+
 end # module
diff --git a/src/Geometry/BoundaryGrids.jl b/src/Geometry/BoundaryGrids.jl
index d8393a40e..58e19bac0 100644
--- a/src/Geometry/BoundaryGrids.jl
+++ b/src/Geometry/BoundaryGrids.jl
@@ -34,8 +34,26 @@ function Triangulation(grid::BoundaryGrid)
   BoundaryTriangulation(trian,grid.descriptor)
 end
 
+function BoundaryGrid(model::DiscreteModel)
+  d = celldim(model)
+  graph = GridGraph(model)
+  facet_to_cells = connections(graph,d-1,d)
+  nfacets = length(facet_to_cells)
+  mask = fill(false,nfacets)
+  _find_on_boundary!(mask,facet_to_cells)
+  BoundaryGrid(model,mask)
+end
+
+function _find_on_boundary!(mask,facet_to_cells)
+  for (facet,cells) in enumerate(facet_to_cells)
+    if length(cells) == 1
+      mask[facet] = true
+    end
+  end
+end
+
 function BoundaryGrid(model::DiscreteModel,tags,icell::Int=1)
-  _tags = _setup_tags(model,tags)
+  _tags = _setup_tags(FaceLabels(model),tags)
   BoundaryGrid(model,_tags,icell)
 end
 
@@ -47,6 +65,12 @@ function BoundaryGrid(model::DiscreteModel,tags::Vector{Int},icell::Int=1)
   nfacets = length(oldfacet_to_label)
   oldfacet_to_mask = fill(false,nfacets)
   _setup_mask!(oldfacet_to_mask,oldfacet_to_label,tag_to_labels,tags)
+  BoundaryGrid(model,oldfacet_to_mask,icell)
+end
+
+function BoundaryGrid(model::DiscreteModel,mask::AbstractVector{Bool},icell::Int=1)
+  d = celldim(model)
+  oldfacet_to_mask = mask
   fgrid = Grid(model,d-1)
   facet_to_oldfacet = findall(oldfacet_to_mask)
   grid = GridPortion(fgrid,facet_to_oldfacet)
@@ -72,6 +96,11 @@ function BoundaryTriangulation(model::DiscreteModel,tags,icell::Int=1)
   Triangulation(grid)
 end
 
+function BoundaryTriangulation(model::DiscreteModel)
+  grid = BoundaryGrid(model)
+  Triangulation(grid)
+end
+
 function _setup_mask!(facet_to_mask,facet_to_label,tag_to_labels,tags)
   nfacets = length(facet_to_mask)
   for facet in 1:nfacets
@@ -96,16 +125,16 @@ function _cell_to_polytope(cell_to_extrussion::ConstantCellValue)
   ConstantCellValue(poly,l)
 end
 
-_setup_tags(model,tags) = tags
+_setup_tags(facelabels,tags) = tags
 
-_setup_tags(model,tag::Integer) = [tag,]
+_setup_tags(facelabels,tag::Integer) = [tag,]
 
-function _setup_tags(model,name::String)
-  _setup_tags(model,[name,])
+function _setup_tags(facelabels,name::String)
+  _setup_tags(facelabels,[name,])
 end
 
-function _setup_tags(model,names::Vector{String})
-  [ tag_from_name(model,s) for s in names ]
+function _setup_tags(facelabels,names::Vector{String})
+  [ tag_from_name(facelabels,s) for s in names ]
 end
 
 end # module
diff --git a/src/Geometry/FaceLabels.jl b/src/Geometry/FaceLabels.jl
index 817063c68..d121584be 100644
--- a/src/Geometry/FaceLabels.jl
+++ b/src/Geometry/FaceLabels.jl
@@ -1,6 +1,7 @@
 module FaceLabelsModule
 
 using Gridap
+using Gridap.Helpers
 using Gridap.CellValuesGallery
 
 export FaceLabels
@@ -10,6 +11,7 @@ export ntags
 export tag_from_name
 export name_from_tag
 export first_tag_on_face
+export add_tag_from_tags!
 
 """
 Classification of nfaces into geometrical and physical labels
@@ -48,6 +50,7 @@ function tag_from_name(l::FaceLabels,name::String)
       return tag
     end
   end
+  @unreachable
   0
 end
 
@@ -73,4 +76,22 @@ function _first_tag_on_face!(face_to_tag,face_to_label,tag_to_labels)
   end
 end
 
+function add_tag_from_tags!(facelabels::FaceLabels, name::String, tags::Vector{Int})
+  @assert ! (name in facelabels.tag_to_name)
+  labels = Int[]
+  for tag in tags
+     for label in labels_on_tag(facelabels,tag)
+       push!(labels,label)
+     end
+  end
+  push!(facelabels.tag_to_labels,sort(collect(Set(labels))))
+  push!(facelabels.tag_to_name,name)
+end
+
+function add_tag_from_tags!(
+  facelabels::FaceLabels, name::String, names::Vector{String})
+  tags = [tag_from_name(facelabels,name) for name in names ]
+  add_tag_from_tags!(facelabels,name,tags)
+end
+
 end #module
diff --git a/src/Geometry/SkeletonGrids.jl b/src/Geometry/SkeletonGrids.jl
index 2bc3e8354..a0c55ae1e 100644
--- a/src/Geometry/SkeletonGrids.jl
+++ b/src/Geometry/SkeletonGrids.jl
@@ -36,7 +36,25 @@ function Triangulation(grid::SkeletonGrid)
   SkeletonTriangulation(trian,grid.descriptor1,grid.descriptor2)
 end
 
-function SkeletonGrid(model::DiscreteModel, tags::Vector{Int})
+function SkeletonGrid(model::DiscreteModel)
+  d = celldim(model)
+  graph = GridGraph(model)
+  facet_to_cells = connections(graph,d-1,d)
+  nfacets = length(facet_to_cells)
+  mask = fill(false,nfacets)
+  _find_on_skeleton!(mask,facet_to_cells)
+  SkeletonGrid(model,mask)
+end
+
+function _find_on_skeleton!(mask,facet_to_cells)
+  for (facet,cells) in enumerate(facet_to_cells)
+    if length(cells) == 2
+      mask[facet] = true
+    end
+  end
+end
+
+function SkeletonGrid(model::DiscreteModel, tags)
   cell1 = 1
   cell2 = 2
   bgrid1 = BoundaryGrid(model,tags,cell1)
@@ -47,9 +65,14 @@ function SkeletonGrid(model::DiscreteModel, tags::Vector{Int})
   SkeletonGrid(grid,descriptor1,descriptor2)
 end
 
-function SkeletonTriangulation(model::DiscreteModel,tags::Vector{Int})
+function SkeletonTriangulation(model::DiscreteModel,tags)
   grid = SkeletonGrid(model,tags)
   Triangulation(grid)
 end
 
+function SkeletonTriangulation(model::DiscreteModel)
+  grid = SkeletonGrid(model)
+  Triangulation(grid)
+end
+
 end # module
diff --git a/src/Integration/NormalVectors.jl b/src/Integration/NormalVectors.jl
index 97e49c02b..e1193d74d 100644
--- a/src/Integration/NormalVectors.jl
+++ b/src/Integration/NormalVectors.jl
@@ -19,6 +19,11 @@ function NormalVector(trian::BoundaryTriangulation)
   NormalVector(desc)
 end
 
+function NormalVector(trian::SkeletonTriangulation)
+  desc = trian.descriptor1
+  NormalVector(desc)
+end
+
 function NormalVector(desc::BoundaryDescriptor)
   phi = desc.cell_phi
   jac = gradient(phi)
diff --git a/src/Integration/SkeletonCellFields.jl b/src/Integration/SkeletonCellFields.jl
index 418da430f..ab24d649f 100644
--- a/src/Integration/SkeletonCellFields.jl
+++ b/src/Integration/SkeletonCellFields.jl
@@ -8,17 +8,25 @@ export SkeletonCellVector
 export SkeletonCellMatrix
 
 export jump
+export mean
 import Gridap: restrict
-import Gridap: gradient
 import Gridap: inner
 import Gridap: integrate
-import Base: -
+import Base: +, -, *
+import Gridap: gradient
+import Gridap: symmetric_gradient
+import Base: div
+import Gridap: trace
+import Gridap: curl
 
 function restrict(
-  cf::IndexCellFieldLike,desc1::BoundaryDescriptor,desc2::BoundaryDescriptor)
+  cf::IndexCellFieldLike{Z,T,N},
+  desc1::BoundaryDescriptor,
+  desc2::BoundaryDescriptor) where {Z,T,N}
+
   cf1 = restrict(cf,desc1)
   cf2 = restrict(cf,desc2)
-  SkeletonPair(cf1,cf2)
+  SkeletonPair{Z-1,T,N}(cf1,cf2)
 end
 
 function restrict(
@@ -28,14 +36,36 @@ function restrict(
 end
 
 struct SkeletonPair{Z,T,N}
-  cellfield1::CellFieldLike{Z,T,N}
-  cellfield2::CellFieldLike{Z,T,N}
+  cellfield1
+  cellfield2
 end
 
-function gradient(sp::SkeletonPair)
-  cf1 = sp.cellfield1
-  cf2 = sp.cellfield2
-  SkeletonPair(gradient(cf1),gradient(cf2))
+for op in (:+,:-,:(gradient),:(symmetric_gradient),:(div),:(trace),:(curl))
+  @eval begin
+    function ($op)(a::SkeletonPair{Z,T,N}) where {Z,T,N}
+      cf1 = a.cellfield1
+      cf2 = a.cellfield2
+      SkeletonPair{Z,T,N}($op(cf1),$op(cf2))
+    end
+  end
+end
+
+for op in (:+, :-, :*)
+  @eval begin
+
+    function ($op)(a::SkeletonPair{Z,T,N},b::CellField) where {Z,T,N}
+      cf1 = a.cellfield1
+      cf2 = a.cellfield2
+      SkeletonPair{Z,T,N}($op(cf1,b),$op(cf2,b))
+    end
+
+    function ($op)(a::CellField,b::SkeletonPair{Z,T,N}) where {Z,T,N}
+      cf1 = b.cellfield1
+      cf2 = b.cellfield2
+      SkeletonPair{Z,T,N}($op(a,cf1),$op(a,cf2))
+    end
+
+  end
 end
 
 function jump(sp::SkeletonPair{Z,T,1}) where {Z,T}
@@ -47,12 +77,24 @@ end
 function jump(sp::SkeletonPair{Z,T,2}) where {Z,T}
   cf1 = sp.cellfield1
   cf2 = sp.cellfield2
-  SkeletonCellBasis(cf1, -cf2)
+  SkeletonCellBasis{Z,T}(cf1, -cf2)
+end
+
+function mean(sp::SkeletonPair{Z,T,1}) where {Z,T}
+  cf1 = sp.cellfield1
+  cf2 = sp.cellfield2
+  0.5*(cf1 + cf2)
+end
+
+function mean(sp::SkeletonPair{Z,T,2}) where {Z,T}
+  cf1 = sp.cellfield1
+  cf2 = sp.cellfield2
+  SkeletonCellBasis{Z,T}(0.5*cf1, 0.5*cf2)
 end
 
 struct SkeletonCellBasis{Z,T}
-  cellbasis1::CellBasis{Z,T}
-  cellbasis2::CellBasis{Z,T}
+  cellbasis1
+  cellbasis2
 end
 
 function inner(a::SkeletonCellBasis{Z},b::CellField{Z}) where Z
@@ -60,7 +102,7 @@ function inner(a::SkeletonCellBasis{Z},b::CellField{Z}) where Z
   cb2 = a.cellbasis2
   cm1 = varinner(cb1,b)
   cm2 = varinner(cb2,b)
-  SkeletonVarinnerVector(cm1,cm2)
+  SkeletonVarinnerVector{Z}(cm1,cm2)
 end
 
 function inner(a::SkeletonCellBasis{Z},b::SkeletonCellBasis{Z}) where Z
@@ -72,43 +114,85 @@ function inner(a::SkeletonCellBasis{Z},b::SkeletonCellBasis{Z}) where Z
   cm12 = varinner(a1,b2)
   cm21 = varinner(a2,b1)
   cm22 = varinner(a2,b2)
-  SkeletonVarinnerMatrix(cm11,cm12,cm21,cm22)
+  SkeletonVarinnerMatrix{Z}(cm11,cm12,cm21,cm22)
 end
 
-struct SkeletonVarinnerVector{D,T}
-  cellmap1::CellMap{Point{D},1,T,2}
-  cellmap2::CellMap{Point{D},1,T,2}
+struct SkeletonVarinnerVector{D}
+  cellmap1
+  cellmap2
 end
 
-function (-)(a::SkeletonVarinnerVector,b::SkeletonVarinnerVector)
-  c1 = a.cellmap1 - b.cellmap1
-  c2 = a.cellmap2 - b.cellmap2
-  SkeletonVarinnerVector(c1,c2)
+for op in (:+, :-)
+  @eval begin
+
+    function ($op)(a::SkeletonVarinnerVector{Z},b::SkeletonVarinnerVector{Z}) where Z
+      c1 = $op(a.cellmap1,b.cellmap1)
+      c2 = $op(a.cellmap2,b.cellmap2)
+      SkeletonVarinnerVector{Z}(c1,c2)
+    end
+    
+    function ($op)(a::SkeletonVarinnerVector{Z}) where Z
+      c1 = $op(a.cellmap1)
+      c2 = $op(a.cellmap2)
+      SkeletonVarinnerVector{Z}(c1,c2)
+    end
+
+  end
 end
 
-function (-)(a::SkeletonVarinnerVector)
-  c1 = - a.cellmap1
-  c2 = - a.cellmap2
-  SkeletonVarinnerVector(c1,c2)
+function (*)(a::Real,b::SkeletonVarinnerVector{Z}) where Z
+  c1 = a*b.cellmap1
+  c2 = a*b.cellmap2
+  SkeletonVarinnerVector{Z}(c1,c2)
 end
 
 struct SkeletonCellVector
-  cellvector1::CellVector
-  cellvector2::CellVector
+  cellvector1
+  cellvector2
+end
+
+struct SkeletonVarinnerMatrix{D}
+  cellmap11
+  cellmap12
+  cellmap21
+  cellmap22
+end
+
+for op in (:+, :-)
+  @eval begin
+
+    function ($op)(a::SkeletonVarinnerMatrix{Z},b::SkeletonVarinnerMatrix{Z}) where Z
+      cellmap11 = $op(a.cellmap11, b.cellmap11)
+      cellmap12 = $op(a.cellmap12, b.cellmap12)
+      cellmap21 = $op(a.cellmap21, b.cellmap21)
+      cellmap22 = $op(a.cellmap22, b.cellmap22)
+      SkeletonVarinnerMatrix{Z}(
+        cellmap11,
+        cellmap12,
+        cellmap21,
+        cellmap22)
+    end
+
+  end
 end
 
-struct SkeletonVarinnerMatrix{D,T}
-  cellmap11::CellMap{Point{D},1,T,3}
-  cellmap12::CellMap{Point{D},1,T,3}
-  cellmap21::CellMap{Point{D},1,T,3}
-  cellmap22::CellMap{Point{D},1,T,3}
+function (*)(a::Real,b::SkeletonVarinnerMatrix{Z}) where Z
+      cellmap11 = a*b.cellmap11
+      cellmap12 = a*b.cellmap12
+      cellmap21 = a*b.cellmap21
+      cellmap22 = a*b.cellmap22
+      SkeletonVarinnerMatrix{Z}(
+        cellmap11,
+        cellmap12,
+        cellmap21,
+        cellmap22)
 end
 
 struct SkeletonCellMatrix
-  cellmatrix11::CellMatrix
-  cellmatrix12::CellMatrix
-  cellmatrix21::CellMatrix
-  cellmatrix22::CellMatrix
+  cellmatrix11
+  cellmatrix12
+  cellmatrix21
+  cellmatrix22
 end
 
 function integrate(
diff --git a/src/Integration/files.jl b/src/Integration/files.jl
index 35c7e5071..fadc9c1fe 100644
--- a/src/Integration/files.jl
+++ b/src/Integration/files.jl
@@ -23,12 +23,12 @@ include("BoundaryDescriptors.jl")
 include("BoundaryTriangulations.jl")
 @reexport using Gridap.BoundaryTriangulations
 
-include("NormalVectors.jl")
-@reexport using Gridap.NormalVectors
-
 include("SkeletonTriangulations.jl")
 @reexport using Gridap.SkeletonTriangulations
 
+include("NormalVectors.jl")
+@reexport using Gridap.NormalVectors
+
 include("BoundaryCellFields.jl")
 @reexport using Gridap.BoundaryCellFields
 
diff --git a/src/MultiField/MultiCellMaps.jl b/src/MultiField/MultiCellMaps.jl
index 1797dec47..075fff819 100644
--- a/src/MultiField/MultiCellMaps.jl
+++ b/src/MultiField/MultiCellMaps.jl
@@ -4,7 +4,7 @@ using Gridap
 
 export MultiCellMap
 import Gridap.CellIntegration: integrate
-import Base: +, -
+import Base: +, -, *
 
 struct MultiCellMap{N}
   blocks::Vector{<:CellMap}
@@ -42,10 +42,20 @@ function (-)(a::MultiCellMap{N},b::MultiCellMap{N}) where N
   MultiCellMap(blocks,fieldids)
 end
 
+function (+)(b::MultiCellMap{N}) where N
+  b
+end
+
 function (-)(b::MultiCellMap{N}) where N
   blocks = CellMap[]
   append!(blocks,[ -k for k in b.blocks])
   MultiCellMap(blocks,b.fieldids)
 end
 
+function (*)(a::Real,b::MultiCellMap{N}) where N
+  blocks = CellMap[]
+  append!(blocks,[ a*k for k in b.blocks])
+  MultiCellMap(blocks,b.fieldids)
+end
+
 end # module
diff --git a/src/MultiField/MultiFEBases.jl b/src/MultiField/MultiFEBases.jl
index 1359d18eb..db8b11c0b 100644
--- a/src/MultiField/MultiFEBases.jl
+++ b/src/MultiField/MultiFEBases.jl
@@ -9,6 +9,7 @@ import Base: div
 import Gridap: trace
 import Gridap: curl
 import Gridap: inner
+import Gridap: varinner
 import Base: +, -, *
 import Base: length, getindex
 import Gridap.FESpaces: FEBasis
@@ -70,6 +71,10 @@ function inner(a::FEBasisWithFieldId,b::CellField)
   MultiCellMap(blocks,fieldids)
 end
 
+function varinner(a::FEBasisWithFieldId,b::CellField)
+  inner(a,b)
+end
+
 function inner(a::FEBasisWithFieldId,f::Function)
   b = CellField(a.febasis.trian,f)
   inner(a,b)
@@ -82,6 +87,21 @@ function inner(a::FEBasisWithFieldId,b::FEBasisWithFieldId)
   MultiCellMap(blocks,fieldids)
 end
 
+function varinner(a::FEBasisWithFieldId,b::FEBasisWithFieldId)
+  inner(a,b)
+end
+
+function restrict(feb::FEBasisWithFieldId,trian::SkeletonTriangulation)
+  sp = restrict(feb.febasis,trian)
+  _new_sp(sp,feb.fieldid)
+end
+
+function _new_sp(sp::SkeletonPair{Z,T,N},fieldid) where {Z,T,N}
+  b1 = FEBasisWithFieldId(sp.cellfield1,fieldid)
+  b2 = FEBasisWithFieldId(sp.cellfield2,fieldid)
+  SkeletonPair{Z,T,N}(b1,b2)
+end
+
 struct MultiFEBasis
   blocks::Vector{<:FEBasisWithFieldId}
 end
@@ -102,9 +122,8 @@ function restrict(mfeb::MultiFEBasis,trian::BoundaryTriangulation)
   MultiFEBasis(blocks)
 end
 
-function restrict(feb::MultiFEBasis,trian::SkeletonTriangulation)
-  @notimplemented
-  # We still need to create a MultiSkeletonPair
+function restrict(mfeb::MultiFEBasis,trian::SkeletonTriangulation)
+  [ restrict(feb,trian) for feb in mfeb.blocks ]
 end
 
 end # module MultiFEBases
diff --git a/src/MultiField/MultiFEFunctions.jl b/src/MultiField/MultiFEFunctions.jl
index e9afd4caf..bc2f327e4 100644
--- a/src/MultiField/MultiFEFunctions.jl
+++ b/src/MultiField/MultiFEFunctions.jl
@@ -52,8 +52,7 @@ function restrict(uh::MultiFEFunction,trian::BoundaryTriangulation)
 end
 
 function restrict(uh::MultiFEFunction,trian::SkeletonTriangulation)
-  @notimplemented
-  # We still need to implement a MultiSkeletonPair
+  [ restrict(ui.cellfield,trian) for ui in uh.fields ]
 end
 
 
diff --git a/test/AlgebraTests/JuliaNLSolversTests.jl b/test/AlgebraTests/JuliaNLSolversTests.jl
new file mode 100644
index 000000000..848d62147
--- /dev/null
+++ b/test/AlgebraTests/JuliaNLSolversTests.jl
@@ -0,0 +1,21 @@
+module JuliaNLSolversTests
+
+using Test
+using Gridap
+
+using ..NonLinearOperatorMocks
+
+op = NonLinearOperatorMock()
+
+nls = JuliaNLSolver(show_trace=false,method=:newton)
+
+x = solve(nls,op)
+@test x ≈ [1.0, 3.0]
+
+cache = solve!(x,nls,op)
+
+solve!(x,nls,op,cache)
+
+cache.result
+
+end # module
diff --git a/test/AlgebraTests/LinearSolversTests.jl b/test/AlgebraTests/LinearSolversTests.jl
index fd8d47a7a..44799db56 100644
--- a/test/AlgebraTests/LinearSolversTests.jl
+++ b/test/AlgebraTests/LinearSolversTests.jl
@@ -28,4 +28,13 @@ b = A*x
 
 test_linear_solver(ls,A,b,x)
 
+n = 10
+A = Laplacian(n,n,1,1)
+x = rand(n^2)
+b = A*x
+
+ls = BackslashSolver()
+
+test_linear_solver(ls,A,b,x)
+
 end # module LinearSolversTests
diff --git a/test/AlgebraTests/runtests.jl b/test/AlgebraTests/runtests.jl
index 95f622486..455d2e4bc 100644
--- a/test/AlgebraTests/runtests.jl
+++ b/test/AlgebraTests/runtests.jl
@@ -8,4 +8,6 @@ include("NonLinearOperatorMocks.jl")
 
 @testset "NonLinearSolvers" begin include("NonLinearSolversTests.jl") end
 
+@testset "JuliaNLSolvers" begin include("JuliaNLSolversTests.jl") end
+
 end # module AlgebraTests
diff --git a/test/FESpacesTests/CLagrangianFESpacesTests.jl b/test/FESpacesTests/CLagrangianFESpacesTests.jl
index c7b0d8c49..86990805a 100644
--- a/test/FESpacesTests/CLagrangianFESpacesTests.jl
+++ b/test/FESpacesTests/CLagrangianFESpacesTests.jl
@@ -160,4 +160,16 @@ fespace = CLagrangianFESpace(T,model,order)
 @test all(fespace.node_and_comp_to_dof .> 0)
 @test fespace.ndiridofs == 0
 
+labels = FaceLabels(model)
+labels = FaceLabels(
+  labels.dim_to_nface_to_label,copy(labels.tag_to_labels),copy(labels.tag_to_name))
+diritag = "diri_boundary"
+add_tag_from_tags!(labels,diritag,diritags)
+dirimask = [true,]
+fespace = CLagrangianFESpace(T,model,labels,order,diritag,dirimask)
+
+
+
+
+
 end # module
diff --git a/test/FESpacesTests/DGFEOperatorsTests.jl b/test/FESpacesTests/DGFEOperatorsTests.jl
new file mode 100644
index 000000000..8b4defea4
--- /dev/null
+++ b/test/FESpacesTests/DGFEOperatorsTests.jl
@@ -0,0 +1,67 @@
+module DGFEOperatorsTests
+
+using Test
+using Gridap
+import Gridap: ∇
+
+u(x) = x[1] + x[2]
+∇u(x) = VectorValue(1.0,1.0)
+∇(::typeof(u)) = ∇u
+f(x) = 0.0
+g(x) = u(x)
+
+L = 1.0
+limits = (0.0, L, 0.0, L)
+ncellx = 4
+model = CartesianDiscreteModel(domain=limits, partition=(ncellx,ncellx))
+
+h = L / ncellx
+
+γ = 10
+
+order = 3
+fespace = DLagrangianFESpace(Float64,model,order)
+
+V = TestFESpace(fespace)
+U = TrialFESpace(fespace)
+
+trian = Triangulation(model)
+quad = CellQuadrature(trian,order=2*order)
+
+btrian = BoundaryTriangulation(model)
+bquad = CellQuadrature(btrian,order=2*order)
+nb = NormalVector(btrian)
+
+strian = SkeletonTriangulation(model)
+squad = CellQuadrature(strian,order=2*order)
+ns = NormalVector(strian)
+
+a_Ω(v,u) = inner(∇(v), ∇(u))
+b_Ω(v) = inner(v,f)
+t_Ω = AffineFETerm(a_Ω,b_Ω,trian,quad)
+
+a_∂Ω(v,u) = (γ/h) * inner(v,u) - inner(v, ∇(u)*nb ) - inner(∇(v)*nb, u)
+b_∂Ω(v) = (γ/h) * inner(v,g) - inner(∇(v)*nb, g)
+t_∂Ω = AffineFETerm(a_∂Ω,b_∂Ω,btrian,bquad)
+
+a_Γ(v,u) = (γ/h) * inner( jump(v*ns), jump(u*ns)) - inner( jump(v*ns), mean(∇(u)) ) - inner( mean(∇(v)), jump(u*ns) ) 
+t_Γ = LinearFETerm(a_Γ,strian,squad)
+
+op = LinearFEOperator(V,U,t_Ω,t_∂Ω,t_Γ)
+
+uh = solve(op)
+
+e = u - uh
+
+#writevtk(trian,"trian",cellfields=["uh"=>uh,"e"=>e])
+
+l2(u) = inner(u,u)
+h1(u) = a_Ω(u,u) + l2(u)
+
+el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
+eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))
+
+@test el2 < 1.e-8
+@test eh1 < 1.e-8
+
+end # module
diff --git a/test/FESpacesTests/DLagrangianFESpacesTests.jl b/test/FESpacesTests/DLagrangianFESpacesTests.jl
index 9d34b4bf0..85b131c16 100644
--- a/test/FESpacesTests/DLagrangianFESpacesTests.jl
+++ b/test/FESpacesTests/DLagrangianFESpacesTests.jl
@@ -127,4 +127,12 @@ fespace = DLagrangianFESpace(T,model,order)
 @test fespace.dirimasks == Tuple{Bool}[]
 @test fespace.ndiridofs == 0
 
+labels = FaceLabels(model)
+labels = FaceLabels(
+  labels.dim_to_nface_to_label,copy(labels.tag_to_labels),copy(labels.tag_to_name))
+diritag = "diri_boundary"
+add_tag_from_tags!(labels,diritag,diritags)
+dirimask = [true,]
+fespace = DLagrangianFESpace(T,model,labels,order,diritag,dirimask)
+
 end # module
diff --git a/test/FESpacesTests/FEBasesTests.jl b/test/FESpacesTests/FEBasesTests.jl
index 864742a2c..a55a063ec 100644
--- a/test/FESpacesTests/FEBasesTests.jl
+++ b/test/FESpacesTests/FEBasesTests.jl
@@ -54,6 +54,17 @@ ca = integrate(cm,btrian,bquad)
 @test isa(ca,CellArray{Float64,1})
 _ = collect(ca)
 
+strian = SkeletonTriangulation(model,"interior")
+squad = CellQuadrature(strian,order=2)
+
+sbh = restrict(bh,strian)
+cm = inner(jump(sbh),mean(sbh))
+ca = integrate(cm,strian,squad)
+
+suh = restrict(uh,strian)
+cm = inner(jump(sbh*mean(suh)),mean(suh))
+ca = integrate(cm,strian,squad)
+
 model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(4,4))
 order = 2
 diritag = "boundary"
diff --git a/test/FESpacesTests/FEOperatorsTestsMixin.jl b/test/FESpacesTests/FEOperatorsTestsMixin.jl
index 52fdc0973..d99306e03 100644
--- a/test/FESpacesTests/FEOperatorsTestsMixin.jl
+++ b/test/FESpacesTests/FEOperatorsTestsMixin.jl
@@ -21,7 +21,7 @@ assem = SparseMatrixAssembler(V,U)
 op = LinearFEOperator(a,b,V,U,assem,trian,quad)
 
 # Define the FESolver
-ls = LUSolver()
+ls = BackslashSolver()
 solver = LinearFESolver(ls)
 
 # Solve!
@@ -70,7 +70,13 @@ gfun(x) = 1.0
 # Construct the FEspace
 order = 1
 diritags = [1,2,3,4,5,6,7]
-fespace = ConformingFESpace(Float64,model,order,diritags)
+labels = FaceLabels(model)
+# This next line is only for testing purposes
+labels = FaceLabels(
+  labels.dim_to_nface_to_label,copy(labels.tag_to_labels),copy(labels.tag_to_name))
+diritag = "diri_boundary"
+add_tag_from_tags!(labels,diritag,diritags)
+fespace = ConformingFESpace(Float64,model,labels,order,diritag)
 
 # Define test and trial
 V = TestFESpace(fespace)
diff --git a/test/FESpacesTests/FESpacesTests.jl b/test/FESpacesTests/FESpacesTests.jl
index 92a216709..0cd2d2e4a 100644
--- a/test/FESpacesTests/FESpacesTests.jl
+++ b/test/FESpacesTests/FESpacesTests.jl
@@ -29,5 +29,13 @@ test_fe_space(V,5,4,mmat,bvec,fun)
 U = TrialFESpace(fespace,fun)
 test_fe_space(U,5,4,mmat,bvec,fun)
 
+fun(x) = x[1] + x[2]
+z = 2.0
+U = TrialFESpace(fespace,z)
+test_fe_space(U,5,4,mmat,bvec,fun)
+
+U = TrialFESpace(fespace,[z,z,z,z])
+test_fe_space(U,5,4,mmat,bvec,fun)
+
 end # module
 
diff --git a/test/FESpacesTests/NonLinearFEOperatorsTests.jl b/test/FESpacesTests/NonLinearFEOperatorsTests.jl
index 8800c8688..0344f0641 100644
--- a/test/FESpacesTests/NonLinearFEOperatorsTests.jl
+++ b/test/FESpacesTests/NonLinearFEOperatorsTests.jl
@@ -46,9 +46,8 @@ op = NonLinearFEOperator(res,jac,V,U,assem,trian,quad)
 
 # Define the FESolver
 ls = LUSolver()
-tol = 1.e-10
-maxiters = 20
-nls = NewtonRaphsonSolver(ls,tol,maxiters)
+nls = JuliaNLSolver(
+  ls,method=:newton,show_trace=false,ftol=1.e-10,iterations=20)
 solver = NonLinearFESolver(nls)
 
 # Solve!
diff --git a/test/FESpacesTests/VectorValuedFEOperatorsTests.jl b/test/FESpacesTests/VectorValuedFEOperatorsTests.jl
index 878e95c76..8eda60d2e 100644
--- a/test/FESpacesTests/VectorValuedFEOperatorsTests.jl
+++ b/test/FESpacesTests/VectorValuedFEOperatorsTests.jl
@@ -3,6 +3,7 @@ module VectorValuedFEOperatorsTests
 using Test
 using Gridap
 import Gridap: ∇
+using LinearAlgebra
 
 # Define manufactured functions
 ufun(x) = VectorValue(x[1] + x[2],x[1])
diff --git a/test/FESpacesTests/runtests.jl b/test/FESpacesTests/runtests.jl
index e15621bf6..a9d4f9b8c 100644
--- a/test/FESpacesTests/runtests.jl
+++ b/test/FESpacesTests/runtests.jl
@@ -24,6 +24,8 @@ using Test
 
 @testset "FEOperators" begin include("FEOperatorsTests.jl") end
 
+@testset "DGFEOperators" begin include("DGFEOperatorsTests.jl") end
+
 @testset "VectorValuedFEOperators" begin include("VectorValuedFEOperatorsTests.jl") end
 
 @testset "NonLinearFEOperators" begin include("NonLinearFEOperatorsTests.jl") end
diff --git a/test/FieldsTests/FieldValuesTests.jl b/test/FieldsTests/FieldValuesTests.jl
index 9b04d073e..7f9290661 100644
--- a/test/FieldsTests/FieldValuesTests.jl
+++ b/test/FieldsTests/FieldValuesTests.jl
@@ -3,6 +3,8 @@ module FieldValuesTests
 using Test
 using Gridap
 using TensorValues
+using LinearAlgebra: dot
+using LinearAlgebra: norm
 
 a = 1
 @test isa(a,FieldValue)
@@ -54,4 +56,9 @@ b = a'
 @test b == TensorValue(1,3,2,4)
 @test a*b == TensorValue(10,14,14,20)
 
+u = VectorValue(1.0,2.0)
+v = VectorValue(2.0,3.0)
+@test dot(u,v) ≈ inner(u,v)
+@test norm(u) ≈ sqrt(inner(u,u))
+
 end # module
diff --git a/test/GeometryTests/FaceLabelsTests.jl b/test/GeometryTests/FaceLabelsTests.jl
index dae37564e..607695d6d 100644
--- a/test/GeometryTests/FaceLabelsTests.jl
+++ b/test/GeometryTests/FaceLabelsTests.jl
@@ -24,4 +24,14 @@ labels = FaceLabels(
 r = [1, 1, 2, 2, 2, 2, 2, 2, 1]
 @test collect(first_tag_on_face(labels,1)) == r
 
+add_tag_from_tags!(labels,"label3",[1,2])
+
+@test labels_on_tag(labels,3) == [1, 2, 3, 4, 5, 6]
+@test tag_from_name(labels,"label3") == 3
+
+add_tag_from_tags!(labels,"label4",["label1","label2"])
+
+@test labels_on_tag(labels,4) == [1, 2, 3, 4, 5, 6]
+@test tag_from_name(labels,"label4") == 4
+
 end # module
diff --git a/test/IntegrationTests/BoundaryTriangulationsTests.jl b/test/IntegrationTests/BoundaryTriangulationsTests.jl
index 931d0372a..0eff2bbd4 100644
--- a/test/IntegrationTests/BoundaryTriangulationsTests.jl
+++ b/test/IntegrationTests/BoundaryTriangulationsTests.jl
@@ -22,4 +22,7 @@ btrian = BoundaryTriangulation(model,tags)
 tags = ["physical_tag_23","physical_tag_24","physical_tag_25"]
 btrian = BoundaryTriangulation(model,tags)
 
+btrian = BoundaryTriangulation(model)
+
+
 end # module
diff --git a/test/IntegrationTests/NormalVectorsTests.jl b/test/IntegrationTests/NormalVectorsTests.jl
index 31bba21fd..80f4cac32 100644
--- a/test/IntegrationTests/NormalVectorsTests.jl
+++ b/test/IntegrationTests/NormalVectorsTests.jl
@@ -26,4 +26,19 @@ x = evaluate(bphi,q)
 #writevtk(btrian,"btrian")
 #writevtk(x,"x",pointdata=["n"=>n_q])
 
+strian = SkeletonTriangulation(model,"interior")
+squad = CellQuadrature(strian,order=2)
+qs = coordinates(squad)
+
+sphi = CellGeomap(strian)
+
+n = NormalVector(strian)
+x = evaluate(sphi,qs)
+n_q = evaluate(n,qs)
+
+#writevtk(strian,"strian")
+#writevtk(x,"x",pointdata=["n"=>n_q])
+
+
+
 end # module
diff --git a/test/IntegrationTests/SkeletonCellFieldsTests.jl b/test/IntegrationTests/SkeletonCellFieldsTests.jl
index 8f508561c..b3a2c7453 100644
--- a/test/IntegrationTests/SkeletonCellFieldsTests.jl
+++ b/test/IntegrationTests/SkeletonCellFieldsTests.jl
@@ -32,18 +32,30 @@ suh = restrict(uh,strian)
 jsu = jump(su)
 jsu_s = evaluate(jsu,s)
 
+msu = mean(su)
+msu_s = evaluate(msu,s)
+
 jsuh = jump(suh)
 jsuh_s = evaluate(jsuh,s)
 
+msuh = mean(suh)
+msuh_s = evaluate(msuh,s)
+
 jsu_grad = jump(∇(su))
 jsu_grad_s = evaluate(jsu_grad,s)
 
+msu_grad = mean(∇(su))
+msu_grad_s = evaluate(msu_grad,s)
+
 jsuh_grad = jump(∇(suh))
 jsuh_grad_s = evaluate(jsuh_grad,s)
 
+msuh_grad = mean(∇(suh))
+msuh_grad_s = evaluate(msuh_grad,s)
+
 #writevtk(strian,"strian",cellfields=[
-# "jumpu"=>jsu,"jumpu_grad"=>jsu_grad,
-# "jumph"=>jsuh,"jumpuh_grad"=>jsuh_grad ])
+# "jumpu"=>jsu,"meanu"=>msu,"jumpu_grad"=>jsu_grad,"meanu_grad"=>msu_grad,
+# "jumph"=>jsuh,"meanh"=>msuh,"jumpuh_grad"=>jsuh_grad,"meanh_grad"=>msuh_grad ])
 
 v = CellBasis(trian)
 du = v
@@ -58,6 +70,12 @@ cv = integrate(cm,strian,squad)
 _ = collect(cv.cellvector1)
 _ = collect(cv.cellvector2)
 
+cm = inner(mean(v_Γ),jump(u_Γ))
+
+cv = integrate(cm,strian,squad)
+_ = collect(cv.cellvector1)
+_ = collect(cv.cellvector2)
+
 cm = inner(jump(v_Γ),jump(du_Γ))
 
 cv = integrate(cm,strian,squad)
@@ -66,6 +84,14 @@ _ = collect(cv.cellmatrix12)
 _ = collect(cv.cellmatrix21)
 _ = collect(cv.cellmatrix22)
 
+cm = inner(mean(v_Γ),jump(du_Γ))
+
+cv = integrate(cm,strian,squad)
+_ = collect(cv.cellmatrix11)
+_ = collect(cv.cellmatrix12)
+_ = collect(cv.cellmatrix21)
+_ = collect(cv.cellmatrix22)
+
 cm = inner(jump(u_Γ),jump(u_Γ))
 cv = integrate(cm,strian,squad)
 _ = collect(cv)
@@ -78,4 +104,8 @@ z = Point(0.0,0.0,0.0)
 zs = fill(z,4)
 @test all([ xi .+ 1 ≈ zs .+ 1 for xi in x ])
 
+x = evaluate(mean(phi_Γ),s)
+x2 = evaluate(phi_Γ.cellfield1,s)
+@test x ≈ x2
+
 end # module
diff --git a/test/IntegrationTests/SkeletonTriangulationsTests.jl b/test/IntegrationTests/SkeletonTriangulationsTests.jl
index 47e2c98c2..a2750709a 100644
--- a/test/IntegrationTests/SkeletonTriangulationsTests.jl
+++ b/test/IntegrationTests/SkeletonTriangulationsTests.jl
@@ -15,4 +15,7 @@ strian = SkeletonTriangulation(trian,descr1,descr2)
 
 test_triangulation(strian)
 
+strian = SkeletonTriangulation(model)
+#writevtk(strian,"strian")
+
 end # module
diff --git a/test/MultiFieldTests/DGMultiFEOperatorsTests.jl b/test/MultiFieldTests/DGMultiFEOperatorsTests.jl
new file mode 100644
index 000000000..14386c077
--- /dev/null
+++ b/test/MultiFieldTests/DGMultiFEOperatorsTests.jl
@@ -0,0 +1,98 @@
+module DGMultiFEOperatorsTests
+
+# Testing Stokes problem with Dirichlet BCs for the velocity
+# and pressure prescribed in a single point
+
+using Test
+using Gridap
+import Gridap: ∇
+
+const T = VectorValue{2,Float64}
+
+# Define manufactured functions
+u1(x) = VectorValue(x[1], x[2])
+u2(x) = x[1] - x[2]
+
+∇u1(x) = TensorValue(1.0,0.0,0.0,1.0)
+
+∇(::typeof(u1)) = ∇u1
+
+b1(x) = VectorValue(1.0,-1.0)
+b2(x) = 2.0
+
+# Construct the discrete model
+model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(2,2))
+
+# Construct the FEspace 1
+order = 2
+diritag = "boundary"
+fespace1 = CLagrangianFESpace(T,model,order,diritag)
+
+# Construct the FEspace 2
+D = 2
+reffe = PDiscRefFE(Float64,D,order-1)
+_fespace2 = DiscFESpace(reffe,model)
+fixeddofs = [1,]
+fespace2 = ConstrainedFESpace(_fespace2,fixeddofs)
+
+# Define test and trial
+V1 = TestFESpace(fespace1)
+V2 = TestFESpace(fespace2)
+V = [V1, V2]
+
+U1 = TrialFESpace(fespace1,u1)
+U2 = TrialFESpace(fespace2)
+U = [U1, U2]
+
+# Define integration mesh and quadrature for volume
+trian = Triangulation(model)
+quad = CellQuadrature(trian,order=2)
+
+strian = SkeletonTriangulation(model)
+squad = CellQuadrature(strian,order=2)
+
+a(v,u) = 
+  inner(∇(v[1]),∇(u[1])) - inner(div(v[1]),u[2]) + inner(v[2],div(u[1]))
+b(v) = inner(v[1],b1) + inner(v[2],b2)
+t_Ω = AffineFETerm(a,b,trian,quad)
+
+# This is a dummy term only for testing purposes. It should not affect the
+# result sine the shape functions for the velocity are continuous (at the moment)
+a_Γ(v,u) = inner(jump(v[1]),jump(u[1]))
+t_Γ = LinearFETerm(a_Γ,strian,squad)
+
+# Define the FEOperator
+op = LinearFEOperator(V,U,t_Ω,t_Γ)
+
+# Solve!
+uh = solve(op)
+
+# Correct the pressure
+A = sum(integrate(u2-uh[2],trian,quad))
+V = sum(integrate((x)->1.0,trian,quad))
+p = uh[2] + A/V
+
+# Define exact solution and error
+e1 = u1 - uh[1]
+
+e2 = u2 - p
+
+#writevtk(trian,"trian",cellfields=["uh2"=>uh[2],"p"=>p])
+
+# Define norms to measure the error
+l2(u) = inner(u,u)
+h1(u) = inner(∇(u),∇(u)) + l2(u)
+
+# Compute errors
+e1l2 = sqrt(sum( integrate(l2(e1),trian,quad) ))
+e1h1 = sqrt(sum( integrate(h1(e1),trian,quad) ))
+
+e2l2 = sqrt(sum( integrate(l2(e2),trian,quad) ))
+
+@test e1l2 < 1.e-8
+@test e1h1 < 1.e-8
+
+@test e2l2 < 1.e-8
+
+
+end # module
diff --git a/test/MultiFieldTests/MultiFEBasesTests.jl b/test/MultiFieldTests/MultiFEBasesTests.jl
index 50743efdb..ecbae848b 100644
--- a/test/MultiFieldTests/MultiFEBasesTests.jl
+++ b/test/MultiFieldTests/MultiFEBasesTests.jl
@@ -104,4 +104,24 @@ ids = ones(Int,ncells(trian))
 cb = CellBasis(trian,σfun,bh[1],ids)
 @test isa(cb,FEBasisWithFieldId)
 
+strian = SkeletonTriangulation(model)
+squad = CellQuadrature(strian,order=2)
+
+sbh = restrict(bh,strian)
+@test isa(sbh[1].cellfield1,FEBasisWithFieldId)
+@test isa(sbh[1].cellfield2,FEBasisWithFieldId)
+
+cm = inner(jump(sbh[1]),jump(sbh[1]))
+ca = integrate(cm,strian,squad)
+
+@test isa(ca.cellmatrix11,MultiCellArray)
+
+uh = zero(V)
+suh = restrict(uh,strian)
+@test isa(suh[1].cellfield1,CellField)
+@test isa(suh[1].cellfield2,CellField)
+
+cm = inner(jump(sbh[1]),jump(suh[1]))
+ca = integrate(cm,strian,squad)
+
 end # module MultiFEBasesTests
diff --git a/test/MultiFieldTests/runtests.jl b/test/MultiFieldTests/runtests.jl
index 9dd2b60d7..732960601 100644
--- a/test/MultiFieldTests/runtests.jl
+++ b/test/MultiFieldTests/runtests.jl
@@ -16,6 +16,8 @@ using Test
 
 @testset "MultiFEOperators" begin include("MultiFEOperatorsTests.jl") end
 
+@testset "DGMultiFEOperators" begin include("DGMultiFEOperatorsTests.jl") end
+
 @testset "VectorValuedMultiFEOperators" begin include("VectorValuedMultiFEOperatorsTests.jl") end
 
 @testset "MultiNonLinearFEOperators" begin include("MultiNonLinearFEOperatorsTests.jl") end