diff --git a/Manifest.toml b/Manifest.toml
index 105e2e9..b92f14d 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -23,11 +23,11 @@ version = "3.3.3"
 [[deps.ArgTools]]
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
 
-[[deps.ArrayInterface]]
-deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
-git-tree-sha1 = "6e8fada11bb015ecf9263f64b156f98b546918c7"
-uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "5.0.5"
+[[deps.ArrayInterfaceCore]]
+deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "d0f59ebfe8d3ea2799fb3fb88742d69978e5843e"
+uuid = "30b0a656-2188-435a-8636-2ec0e6a096e2"
+version = "0.1.10"
 
 [[deps.ArrayLayouts]]
 deps = ["FillArrays", "LinearAlgebra", "SparseArrays"]
@@ -72,15 +72,15 @@ version = "1.16.1+1"
 
 [[deps.ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "9950387274246d08af38f6eef8cb5480862a435f"
+git-tree-sha1 = "9489214b993cd42d17f44c36e359bf6a7c919abf"
 uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "1.14.0"
+version = "1.15.0"
 
 [[deps.ChangesOfVariables]]
 deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
-git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
+git-tree-sha1 = "1e315e3f4b0b7ce40feded39c73049692126cf53"
 uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
-version = "0.1.2"
+version = "0.1.3"
 
 [[deps.CodecZlib]]
 deps = ["TranscodingStreams", "Zlib_jll"]
@@ -89,16 +89,22 @@ uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
 version = "0.7.0"
 
 [[deps.ColorSchemes]]
-deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
-git-tree-sha1 = "12fc73e5e0af68ad3137b886e3f7c1eacfca2640"
+deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Random"]
+git-tree-sha1 = "7297381ccb5df764549818d9a7d57e45f1057d30"
 uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
-version = "3.17.1"
+version = "3.18.0"
 
 [[deps.ColorTypes]]
 deps = ["FixedPointNumbers", "Random"]
-git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
+git-tree-sha1 = "0f4e115f6f34bbe43c19751c90a38b2f380637b9"
 uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
-version = "0.11.0"
+version = "0.11.3"
+
+[[deps.ColorVectorSpace]]
+deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"]
+git-tree-sha1 = "d08c20eef1f2cbc6e60fd3612ac4340b89fea322"
+uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
+version = "0.9.9"
 
 [[deps.Colors]]
 deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
@@ -119,9 +125,9 @@ version = "0.3.0"
 
 [[deps.Compat]]
 deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
-git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
+git-tree-sha1 = "9be8be1d8a6f44b96482c8af52238ea7987da3e3"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "3.42.0"
+version = "3.45.0"
 
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -139,21 +145,21 @@ uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
 version = "4.1.1"
 
 [[deps.DataAPI]]
-git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
+git-tree-sha1 = "fb5f5316dd3fd4c5e7c30a24d50643b73e37cd40"
 uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
-version = "1.9.0"
+version = "1.10.0"
 
 [[deps.DataFrames]]
 deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
-git-tree-sha1 = "ae02104e835f219b8930c7664b8012c93475c340"
+git-tree-sha1 = "daa21eb85147f72e41f6352a57fccea377e310a9"
 uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-version = "1.3.2"
+version = "1.3.4"
 
 [[deps.DataStructures]]
 deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
-git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
+git-tree-sha1 = "d1fff3a548102f48987a52a2e0d114fa97d730f0"
 uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-version = "0.18.11"
+version = "0.18.13"
 
 [[deps.DataValueInterfaces]]
 git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
@@ -182,9 +188,9 @@ version = "1.0.3"
 
 [[deps.DiffRules]]
 deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
-git-tree-sha1 = "dd933c4ef7b4c270aacd4eb88fa64c147492acf0"
+git-tree-sha1 = "28d605d9a0ac17118fe2c5e9ce0fbb76c3ceb120"
 uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
-version = "1.10.0"
+version = "1.11.0"
 
 [[deps.Distances]]
 deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"]
@@ -218,16 +224,11 @@ git-tree-sha1 = "3f3a2501fa7236e9b911e0f7a588c657e822bb6d"
 uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
 version = "2.2.3+0"
 
-[[deps.EllipsisNotation]]
-git-tree-sha1 = "18ee049accec8763be17a933737c1dd0fdf8673a"
-uuid = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
-version = "1.0.0"
-
 [[deps.Expat_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "ae13fcbc7ab8f16b0856729b050ef0c446aa3492"
+git-tree-sha1 = "bad72f730e9e91c08d9427d5e8db95478a3c323d"
 uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
-version = "2.4.4+0"
+version = "2.4.8+0"
 
 [[deps.FFMPEG]]
 deps = ["FFMPEG_jll"]
@@ -249,9 +250,9 @@ version = "0.4.9"
 
 [[deps.FileIO]]
 deps = ["Pkg", "Requires", "UUIDs"]
-git-tree-sha1 = "80ced645013a5dbdc52cf70329399c35ce007fae"
+git-tree-sha1 = "9267e5f50b0e12fdfd5a2455534345c4cf2c7f7a"
 uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.13.0"
+version = "1.14.0"
 
 [[deps.FillArrays]]
 deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
@@ -260,10 +261,10 @@ uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
 version = "0.12.8"
 
 [[deps.FiniteDiff]]
-deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"]
-git-tree-sha1 = "56956d1e4c1221000b7781104c58c34019792951"
+deps = ["ArrayInterfaceCore", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"]
+git-tree-sha1 = "a0700c21266b55bf62c22e75af5668aa7841b500"
 uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
-version = "2.11.0"
+version = "2.12.1"
 
 [[deps.FixedPointNumbers]]
 deps = ["Statistics"]
@@ -291,9 +292,9 @@ version = "0.4.2"
 
 [[deps.ForwardDiff]]
 deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
-git-tree-sha1 = "1bd6fc0c344fc0cbee1f42f8d2e7ec8253dda2d2"
+git-tree-sha1 = "2f18915445b248731ec5db4e4a17e451020bf21e"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "0.10.25"
+version = "0.10.30"
 
 [[deps.FreeType2_jll]]
 deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
@@ -319,15 +320,15 @@ version = "3.3.6+0"
 
 [[deps.GR]]
 deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "RelocatableFolders", "Serialization", "Sockets", "Test", "UUIDs"]
-git-tree-sha1 = "9f836fb62492f4b0f0d3b06f55983f2704ed0883"
+git-tree-sha1 = "c98aea696662d09e215ef7cda5296024a9646c75"
 uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
-version = "0.64.0"
+version = "0.64.4"
 
 [[deps.GR_jll]]
 deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"]
-git-tree-sha1 = "a6c850d77ad5118ad3be4bd188919ce97fffac47"
+git-tree-sha1 = "3a233eeeb2ca45842fe100e0413936834215abf5"
 uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
-version = "0.64.0+0"
+version = "0.64.4+0"
 
 [[deps.GTK3_jll]]
 deps = ["ATK_jll", "Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "FriBidi_jll", "Glib_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Libepoxy_jll", "Pango_jll", "Pkg", "Wayland_jll", "Xorg_libX11_jll", "Xorg_libXcomposite_jll", "Xorg_libXcursor_jll", "Xorg_libXdamage_jll", "Xorg_libXext_jll", "Xorg_libXfixes_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll", "Xorg_libXrender_jll", "at_spi2_atk_jll", "gdk_pixbuf_jll", "iso_codes_jll", "xkbcommon_jll"]
@@ -355,9 +356,9 @@ version = "2.68.3+2"
 
 [[deps.Graphics]]
 deps = ["Colors", "LinearAlgebra", "NaNMath"]
-git-tree-sha1 = "1c5a84319923bea76fa145d49e93aa4394c73fc2"
+git-tree-sha1 = "d61890399bc535850c4bf08e4e0d3a7ad0f21cbd"
 uuid = "a2bd30eb-e257-5431-a919-1863eab51364"
-version = "1.1.1"
+version = "1.1.2"
 
 [[deps.Graphite2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -373,6 +374,12 @@ repo-url = "https://github.com/gridap/Gridap.jl"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 version = "0.17.12"
 
+[[deps.GridapDistributed]]
+deps = ["FillArrays", "Gridap", "LinearAlgebra", "MPI", "PartitionedArrays", "SparseArrays", "SparseMatricesCSR", "WriteVTK"]
+git-tree-sha1 = "0b67924b214c9735db46071646e6ab1e1817eddd"
+uuid = "f9701e48-63b3-45aa-9a63-9bc6c271f355"
+version = "0.2.6"
+
 [[deps.Grisu]]
 git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
 uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
@@ -386,9 +393,9 @@ version = "1.2.1"
 
 [[deps.GtkObservables]]
 deps = ["Cairo", "Colors", "Dates", "FixedPointNumbers", "Graphics", "Gtk", "IntervalSets", "LinearAlgebra", "Observables", "Reexport", "RoundingIntegers"]
-git-tree-sha1 = "96dcca8d49566c6a0fddc98fb46f7b810de6ef63"
+git-tree-sha1 = "cf87f031fee932b90023ea37207c7a1de8caee6f"
 uuid = "8710efd8-4ad6-11eb-33ea-2d5ceb25a41c"
-version = "1.2.0"
+version = "1.2.3"
 
 [[deps.HTTP]]
 deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
@@ -402,11 +409,6 @@ git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3"
 uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
 version = "2.8.1+1"
 
-[[deps.IfElse]]
-git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
-uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
-version = "0.1.1"
-
 [[deps.IndirectArrays]]
 git-tree-sha1 = "012e604e1c7458645cb8b436f8fba789a51b257f"
 uuid = "9b13fd28-a010-5f03-acff-a1bbcff69959"
@@ -422,16 +424,16 @@ deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
 [[deps.IntervalSets]]
-deps = ["Dates", "EllipsisNotation", "Statistics"]
-git-tree-sha1 = "bcf640979ee55b652f3b01650444eb7bbe3ea837"
+deps = ["Dates", "Statistics"]
+git-tree-sha1 = "ad841eddfb05f6d9be0bff1fa48dcae32f134a2d"
 uuid = "8197267c-284f-5f27-9208-e0e47529a953"
-version = "0.5.4"
+version = "0.6.2"
 
 [[deps.InverseFunctions]]
 deps = ["Test"]
-git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
+git-tree-sha1 = "c6cf981474e7094ce044168d329274d797843467"
 uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
-version = "0.1.3"
+version = "0.1.6"
 
 [[deps.InvertedIndices]]
 git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
@@ -448,6 +450,12 @@ git-tree-sha1 = "fa6287a4469f5e048d763df38279ee729fbd44e5"
 uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 version = "1.4.0"
 
+[[deps.IterativeSolvers]]
+deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays"]
+git-tree-sha1 = "1169632f425f79429f245113b775a0e3d121457c"
+uuid = "42fd0dbc-a981-5370-80f2-aaf504508153"
+version = "0.9.2"
+
 [[deps.IteratorInterfaceExtensions]]
 git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
 uuid = "82899510-4779-5014-852e-03e436cf321d"
@@ -502,9 +510,13 @@ version = "1.3.0"
 
 [[deps.Latexify]]
 deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
-git-tree-sha1 = "4f00cc36fede3c04b8acf9b2e2763decfdcecfa6"
+git-tree-sha1 = "46a39b9c58749eefb5f2dc1178cb8fab5332b1ab"
 uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
-version = "0.15.13"
+version = "0.15.15"
+
+[[deps.LazyArtifacts]]
+deps = ["Artifacts", "Pkg"]
+uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 
 [[deps.LeftChildRightSiblingTrees]]
 deps = ["AbstractTrees"]
@@ -581,9 +593,9 @@ version = "2.52.4+0"
 
 [[deps.Libtiff_jll]]
 deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
-git-tree-sha1 = "c9551dd26e31ab17b86cbd00c2ede019c08758eb"
+git-tree-sha1 = "3eb79b0ca5764d4799c06699573fd8f533259713"
 uuid = "89763e89-9b03-5906-acba-b20f662cd828"
-version = "4.3.0+1"
+version = "4.4.0+0"
 
 [[deps.Libuuid_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -609,13 +621,25 @@ uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [[deps.LogExpFunctions]]
 deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
-git-tree-sha1 = "58f25e56b706f95125dcb796f39e1fb01d913a71"
+git-tree-sha1 = "09e4b894ce6a976c354a69041a04748180d43637"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-version = "0.3.10"
+version = "0.3.15"
 
 [[deps.Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
+[[deps.MPI]]
+deps = ["Distributed", "DocStringExtensions", "Libdl", "MPICH_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "Pkg", "Random", "Requires", "Serialization", "Sockets"]
+git-tree-sha1 = "d56a80d8cf8b9dc3050116346b3d83432b1912c0"
+uuid = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+version = "0.19.2"
+
+[[deps.MPICH_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "9f781ffc4020c24ca05d83fa4eb4df0d34832b18"
+uuid = "7cb0a576-ebde-5e09-9194-50597f1243b4"
+version = "4.0.2+1"
+
 [[deps.MacroTools]]
 deps = ["Markdown", "Random"]
 git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
@@ -647,6 +671,12 @@ git-tree-sha1 = "81e123ea81d6081fe8b733dbe79e1291d55cfb0f"
 uuid = "85b6ec6f-f7df-4429-9514-a64bcd9ee824"
 version = "0.4.6"
 
+[[deps.MicrosoftMPI_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "a16aa086d335ed7e0170c5265247db29172af2f9"
+uuid = "9237b28f-5490-5468-be7b-bb81f5f5e6cf"
+version = "10.1.3+2"
+
 [[deps.Missings]]
 deps = ["DataAPI"]
 git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
@@ -686,9 +716,9 @@ version = "0.4.10"
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 
 [[deps.Observables]]
-git-tree-sha1 = "fe29afdef3d0c4a8286128d4e45cc50621b1e43d"
+git-tree-sha1 = "dfd8d34871bc3ad08cd16026c1828e271d554db9"
 uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
-version = "0.4.0"
+version = "0.5.1"
 
 [[deps.Ogg_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -704,6 +734,12 @@ uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
 deps = ["Artifacts", "Libdl"]
 uuid = "05823500-19ac-5b8b-9628-191a04bc5112"
 
+[[deps.OpenMPI_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
+git-tree-sha1 = "19ec7d0311aa5fb5fe537dc6eeaec86942b64caf"
+uuid = "fe0851c0-eecd-5654-98d4-656369965a5c"
+version = "4.1.3+0"
+
 [[deps.OpenSSL_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "ab05aa4cc89736e95915b01e7279e61b1bfe33b8"
@@ -747,9 +783,15 @@ version = "0.12.3"
 
 [[deps.Parsers]]
 deps = ["Dates"]
-git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
+git-tree-sha1 = "1285416549ccfcdf0c50d4997a94331e88d68413"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.2.3"
+version = "2.3.1"
+
+[[deps.PartitionedArrays]]
+deps = ["Distances", "IterativeSolvers", "LinearAlgebra", "MPI", "Printf", "SparseArrays", "SparseMatricesCSR"]
+git-tree-sha1 = "88ff2293fd57089a4036a3056ba058ae9806111b"
+uuid = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
+version = "0.2.10"
 
 [[deps.Pixman_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -762,10 +804,10 @@ deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markd
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
 [[deps.PlotThemes]]
-deps = ["PlotUtils", "Requires", "Statistics"]
-git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
+deps = ["PlotUtils", "Statistics"]
+git-tree-sha1 = "8162b2f8547bc23876edd0c5181b27702ae58dce"
 uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
-version = "2.0.1"
+version = "3.0.0"
 
 [[deps.PlotUtils]]
 deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
@@ -775,21 +817,21 @@ version = "1.2.0"
 
 [[deps.Plots]]
 deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "Pkg", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun", "Unzip"]
-git-tree-sha1 = "5f6e1309595e95db24342e56cd4dabd2159e0b79"
+git-tree-sha1 = "9e42de869561d6bdf8602c57ec557d43538a92f0"
 uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-version = "1.27.3"
+version = "1.29.1"
 
 [[deps.PooledArrays]]
 deps = ["DataAPI", "Future"]
-git-tree-sha1 = "28ef6c7ce353f0b35d0df0d5930e0d072c1f5b9b"
+git-tree-sha1 = "a6062fe4063cdafe78f4a0a81cfffb89721b30e7"
 uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
-version = "1.4.1"
+version = "1.4.2"
 
 [[deps.Preferences]]
 deps = ["TOML"]
-git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
+git-tree-sha1 = "47e5f437cc0e7ef2ce8406ce1e7e24d44915f88d"
 uuid = "21216c6a-2e73-6563-6e65-726566657250"
-version = "1.2.5"
+version = "1.3.0"
 
 [[deps.PrettyTables]]
 deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
@@ -807,15 +849,15 @@ uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
 
 [[deps.ProfileView]]
 deps = ["Cairo", "Colors", "FileIO", "FlameGraphs", "Graphics", "Gtk", "GtkObservables", "InteractiveUtils", "IntervalSets", "MethodAnalysis", "Preferences", "Profile", "UUIDs"]
-git-tree-sha1 = "1638391c64d5ae2853b54823df3f98150dc202f3"
+git-tree-sha1 = "17c95da7223ca01bccedde748ae4ccfebb494106"
 uuid = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7"
-version = "1.5.0"
+version = "1.5.1"
 
 [[deps.Qt5Base_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"]
-git-tree-sha1 = "ad368663a5e20dbb8d6dc2fddeefe4dae0781ae8"
+git-tree-sha1 = "c6c0f690d0cc7caddb74cef7aa847b824a16b256"
 uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1"
-version = "5.15.3+0"
+version = "5.15.3+1"
 
 [[deps.QuadGK]]
 deps = ["DataStructures", "LinearAlgebra"]
@@ -907,21 +949,15 @@ version = "0.6.6"
 
 [[deps.SpecialFunctions]]
 deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
+git-tree-sha1 = "a9e798cae4867e3a41cae2dd9eb60c047f1212db"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.1.4"
-
-[[deps.Static]]
-deps = ["IfElse"]
-git-tree-sha1 = "87e9954dfa33fd145694e42337bdd3d5b07021a6"
-uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
-version = "0.6.0"
+version = "2.1.6"
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "Random", "Statistics"]
-git-tree-sha1 = "4f6ec5d99a28e1a749559ef7dd518663c5eca3d5"
+git-tree-sha1 = "383a578bdf6e6721f480e749d503ebc8405a0b22"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.4.3"
+version = "1.4.6"
 
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
@@ -929,9 +965,9 @@ uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[deps.StatsAPI]]
 deps = ["LinearAlgebra"]
-git-tree-sha1 = "c3d8ba7f3fa0625b062b82853a7d5229cb728b6b"
+git-tree-sha1 = "2c11d7290036fe7aac9038ff312d3b3a2a5bf89e"
 uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
-version = "1.2.1"
+version = "1.4.0"
 
 [[deps.StatsBase]]
 deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
@@ -941,9 +977,9 @@ version = "0.33.16"
 
 [[deps.StructArrays]]
 deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
-git-tree-sha1 = "57617b34fa34f91d536eb265df67c2d4519b8b98"
+git-tree-sha1 = "9abba8f8fb8458e9adf07c8a2377a070674a24f1"
 uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
-version = "0.6.5"
+version = "0.6.8"
 
 [[deps.SuiteSparse]]
 deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
@@ -969,6 +1005,12 @@ version = "1.7.0"
 deps = ["ArgTools", "SHA"]
 uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
 
+[[deps.TensorCore]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6"
+uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
+version = "0.1.1"
+
 [[deps.Test]]
 deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
@@ -1027,9 +1069,9 @@ version = "1.14.2"
 
 [[deps.XML2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
+git-tree-sha1 = "58443b63fb7e465a8a7210828c91c08b92132dff"
 uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
-version = "2.9.12+0"
+version = "2.9.14+0"
 
 [[deps.XSLT_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
diff --git a/Project.toml b/Project.toml
index 3326394..cd1d468 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,12 +8,16 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DrWatson = "634d3b9d-ee7a-5ddf-bec9-22491ea816e1"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
+GridapDistributed = "f9701e48-63b3-45aa-9a63-9bc6c271f355"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ProfileView = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
+GridapDistributed = "0.2.6"
 DrWatson = "1"
 Gridap = "0.17"
diff --git a/src/Distributed/HybridAffineFEOperators.jl b/src/Distributed/HybridAffineFEOperators.jl
new file mode 100644
index 0000000..4a7282f
--- /dev/null
+++ b/src/Distributed/HybridAffineFEOperators.jl
@@ -0,0 +1,85 @@
+function _merge_bulk_and_skeleton_contributions(
+    matcontribs::GridapDistributed.DistributedDomainContribution,
+    veccontribs::GridapDistributed.DistributedDomainContribution)
+  dmat,dvec=map_parts(matcontribs.contribs,veccontribs.contribs) do matseq, vecseq
+    _merge_bulk_and_skeleton_contributions(matseq,vecseq)
+  end
+  GridapDistributed.DistributedDomainContribution(dmat),
+          GridapDistributed.DistributedDomainContribution(dvec)
+end
+
+function Gridap.FESpaces._pair_contribution_when_possible(
+    obiform::GridapDistributed.DistributedDomainContribution,
+    oliform::GridapDistributed.DistributedDomainContribution)
+  dmv,dm,dv=map_parts(obiform.contribs,oliform.contribs) do obiform, oliform
+    Gridap.FESpaces._pair_contribution_when_possible(obiform,oliform)
+  end
+  GridapDistributed.DistributedDomainContribution(dmv),
+    GridapDistributed.DistributedDomainContribution(dm),
+      GridapDistributed.DistributedDomainContribution(dv)
+end
+
+function Gridap.FESpaces._collect_cell_matrix_and_vector(
+  trial  :: GridapDistributed.DistributedFESpace,
+  test   :: GridapDistributed.DistributedFESpace,
+  matvec :: GridapDistributed.DistributedDomainContribution,
+  mat    :: GridapDistributed.DistributedDomainContribution,
+  vec    :: GridapDistributed.DistributedDomainContribution)
+  map_parts(GridapDistributed.local_views(trial),
+            GridapDistributed.local_views(test),
+            matvec.contribs,mat.contribs,vec.contribs) do trial, test, matvec, mat, vec
+     Gridap.FESpaces._collect_cell_matrix_and_vector(trial,test,matvec,mat,vec)
+  end
+end
+
+function _add_static_condensation(
+   matvec::GridapDistributed.DistributedDomainContribution,
+   bulk_fields, skeleton_fields)
+   dmv=map_parts(matvec.contribs) do matvec
+    _add_static_condensation(matvec,bulk_fields,skeleton_fields)
+  end
+  GridapDistributed.DistributedDomainContribution(dmv)
+end
+
+function Gridap.FESpaces._attach_dirichlet(
+  matvec:: GridapDistributed.DistributedDomainContribution,
+  mat   :: GridapDistributed.DistributedDomainContribution,
+  uhd)
+  dmv,dm=map_parts(matvec.contribs,mat.contribs,uhd.fields) do matvec, mat, uhd
+    m,v=Gridap.FESpaces._attach_dirichlet(matvec,mat, uhd)
+    m,v
+  end
+  GridapDistributed.DistributedDomainContribution(dmv),
+      GridapDistributed.DistributedDomainContribution(dm)
+end
+
+
+function _compute_hybridizable_from_skeleton_free_dof_values(
+    skeleton_fe_function :: GridapDistributed.DistributedCellField,
+    trial_hybridizable   :: GridapDistributed.DistributedFESpace,
+    test_hybridizable    :: GridapDistributed.DistributedFESpace,
+    trial_skeleton       :: GridapDistributed.DistributedFESpace,
+    matvec               :: GridapDistributed.DistributedDomainContribution,
+    bulk_fields, skeleton_fields)
+
+    function f(skel_fe, trial_hyb, test_hyb, trial_skel, matvec)
+      _compute_hybridizable_from_skeleton_free_dof_values(skel_fe,
+                                                          trial_hyb,
+                                                          test_hyb,
+                                                          trial_skel,
+                                                          matvec,
+                                                          bulk_fields,
+                                                          skeleton_fields)
+    end
+    # f(skeleton_fe_function.fields.parts[1],
+    #   trial_hybridizable.part_fe_space.parts[1],
+    #   test_hybridizable.part_fe_space.parts[1],
+    #   trial_skeleton.spaces.parts[1],
+    #   matvec.contribs.parts[1])
+    values=map_parts(f,GridapDistributed.local_views(skeleton_fe_function),
+             GridapDistributed.local_views(trial_hybridizable),
+             GridapDistributed.local_views(test_hybridizable),
+             GridapDistributed.local_views(trial_skeleton),
+             GridapDistributed.local_views(matvec))
+    PVector(values,trial_hybridizable.gids)
+end
diff --git a/src/Distributed/Skeleton.jl b/src/Distributed/Skeleton.jl
new file mode 100644
index 0000000..e46b493
--- /dev/null
+++ b/src/Distributed/Skeleton.jl
@@ -0,0 +1,49 @@
+
+function SkeletonTriangulation(model::GridapDistributed.DistributedDiscreteModel)
+  SkeletonTriangulation(no_ghost,model::GridapDistributed.DistributedDiscreteModel)
+end
+
+
+function _sign_flips(model::GridapDistributed.DistributedDiscreteModel)
+  cell_reffes = map_parts(model.models) do m
+    basis,reffe_args,reffe_kwargs = ReferenceFE(raviart_thomas,Float64,0)
+    reffe = ReferenceFE(m,basis,reffe_args...;reffe_kwargs...)
+  end
+  GridapDistributed._generate_sign_flips(model,cell_reffes)
+end
+
+function SkeletonTriangulation(portion::GridapDistributed.WithGhost,
+                               model::GridapDistributed.DistributedDiscreteModel)
+   dtrians=map_parts(model.models,_sign_flips(model)) do model, sign_flip
+      Skeleton(model,sign_flip)
+   end
+   GridapDistributed.DistributedTriangulation(dtrians,model)
+end
+
+function SkeletonTriangulation(portion::GridapDistributed.NoGhost,
+                               model::GridapDistributed.DistributedDiscreteModel)
+  cell_gids=get_cell_gids(model)
+  dtrians=map_parts(model.models,
+                    _sign_flips(model),
+                    cell_gids.partition) do model, sign_flip, partition
+    cell_to_parent_cell=
+      findall([partition.lid_to_part[cell]==partition.part
+                 for cell=1:length(partition.lid_to_part)])
+    Skeleton(model,cell_to_parent_cell,sign_flip)
+  end
+  GridapDistributed.DistributedTriangulation(dtrians,model)
+end
+
+function get_cell_normal_vector(dskeleton::GridapDistributed.DistributedTriangulation)
+  fields=map_parts(dskeleton.trians) do skeleton
+    get_cell_normal_vector(skeleton)
+  end
+  GridapDistributed.DistributedCellField(fields)
+end
+
+function get_cell_owner_normal_vector(dskeleton::GridapDistributed.DistributedTriangulation)
+  fields=map_parts(dskeleton.trians) do skeleton
+    get_cell_owner_normal_vector(skeleton)
+  end
+  GridapDistributed.DistributedCellField(fields)
+end
diff --git a/src/GridapHybrid.jl b/src/GridapHybrid.jl
index 3f0166d..db00784 100644
--- a/src/GridapHybrid.jl
+++ b/src/GridapHybrid.jl
@@ -37,4 +37,9 @@ include("HybridLinearSolvers.jl")
 include("GridapAPIExtensions.jl")
 include("GridapTmpModifications.jl")
 
+using GridapDistributed
+using PartitionedArrays
+include("Distributed/Skeleton.jl")
+include("Distributed/HybridAffineFEOperators.jl")
+
 end # module
diff --git a/src/GridapTmpModifications.jl b/src/GridapTmpModifications.jl
index 837231a..8b2a3c8 100644
--- a/src/GridapTmpModifications.jl
+++ b/src/GridapTmpModifications.jl
@@ -302,11 +302,12 @@ function Gridap.Arrays.evaluate!(cache,
 end
 
 function _generate_cell_lface_dofs_from_cell_dofs(glue,
+                                                  cell_grid,
                                                   cell_wise_facets,
                                                   cell_values_field,
                                                   facets_dofs_ids)
   facet_sizes=lazy_map(x->length(x),facets_dofs_ids)
-  cell_wise_facet_sizes=SkeletonVectorFromFacetVector(glue,cell_wise_facets,facet_sizes)
+  cell_wise_facet_sizes=SkeletonVectorFromFacetVector(glue,cell_grid,cell_wise_facets,facet_sizes)
   m=SkeletonVectorFromNonBlockedSkeletonVector(cell_wise_facet_sizes)
   lazy_map(m,cell_values_field,collect(1:length(cell_values_field)))
 end
@@ -335,13 +336,15 @@ function Gridap.FESpaces._change_argument(
         cell_lface_dof_values =
          _generate_cell_lface_dofs_from_cell_dofs(
             trian.glue,
+            trian.grid.parent,
             cell_wise_facets,
             cell_values_field,
             get_cell_dof_ids(Ui))
         Ui_basis = get_fe_basis(Ui)
-        Ui_basis_data = get_data(Ui_basis)
+        Ui_basis_data = Gridap.CellData.get_data(Ui_basis)
         Ui_basis_cell_lface_data = SkeletonVectorFromFacetVector(
           trian.glue,
+          trian.grid.parent,
           cell_wise_facets,
           Ui_basis_data)
         cell_field = lazy_map(linear_combination, cell_lface_dof_values, Ui_basis_cell_lface_data)
diff --git a/src/HybridAffineFEOperators.jl b/src/HybridAffineFEOperators.jl
index 3edf305..64ba799 100644
--- a/src/HybridAffineFEOperators.jl
+++ b/src/HybridAffineFEOperators.jl
@@ -1,7 +1,7 @@
 struct HybridAffineFEOperator{TB,TS} <: FEOperator
     weakform::Function
-    trial::MultiFieldFESpace
-    test::MultiFieldFESpace
+    trial
+    test
     bulk_fields::TB
     skeleton_fields::TS
     skeleton_op::AffineFEOperator
@@ -9,8 +9,8 @@ end
 
 function HybridAffineFEOperator(
   weakform::Function,
-  trial::MultiFieldFESpace,
-  test::MultiFieldFESpace,
+  trial,
+  test,
   bulk_fields::TB,
   skeleton_fields::TS) where {TB <: Vector{<:Integer},TS <: Vector{<:Integer}}
 
@@ -27,6 +27,8 @@ function HybridAffineFEOperator(
   # Pair LHS and RHS terms associated to SkeletonTriangulation
     matvec, mat, vec = Gridap.FESpaces._pair_contribution_when_possible(obiform, oliform)
 
+
+
   # Add StaticCondensationMap to matvec terms
     matvec = _add_static_condensation(matvec, bulk_fields, skeleton_fields)
 
@@ -71,7 +73,7 @@ end
 
 function Gridap.FESpaces.solve!(uh, solver::LinearFESolver, op::HybridAffineFEOperator, cache)
   # Solve linear system defined on the skeleton
-    lh = solve(op.skeleton_op)
+    lh = solve(solver,op.skeleton_op)
 
   # Invoke weak form of the hybridizable system
     u = get_trial_fe_basis(op.trial)
@@ -121,25 +123,39 @@ function _compute_hybridizable_from_skeleton_free_dof_values(skeleton_fe_functio
     cell_wise_facets = _get_cell_wise_facets(model)
     cells_around_facets = _get_cells_around_facets(model)
 
+    if (isa(Γ.grid.parent,Gridap.Geometry.GridView))
+      cell_to_parent_cell = Γ.grid.parent.cell_to_parent_cell
+      ifacet_to_facet, facet_to_icell =
+       _find_faces_touched_by_cells(cell_to_parent_cell,
+          cell_wise_facets,
+          cells_around_facets)
+    end
+
     nfields = length(bulk_fields) + length(skeleton_fields)
     m = Gridap.Fields.BlockMap(nfields, skeleton_fields)
     L = Gridap.FESpaces.get_fe_space(skeleton_fe_function)
+    L_cell_dof_ids = get_cell_dof_ids(L)
+    if (isa(Γ.grid.parent,Gridap.Geometry.GridView))
+      L_cell_dof_ids=lazy_map(Reindex(L_cell_dof_ids),ifacet_to_facet)
+    end
+
     lhₑ = lazy_map(m,
                  convert_cell_wise_dofs_array_to_facet_dofs_array(
                  cells_around_facets,
                  cell_wise_facets,
                  lhₖ,
-                 get_cell_dof_ids(L))...)
+                 L_cell_dof_ids)...)
 
     assem = SparseMatrixAssembler(trial_hybridizable, test_hybridizable)
-    lhₑ_dofs = get_cell_dof_ids(trial_hybridizable, get_triangulation(L))
+    Ltrian=get_triangulation(L)
+    if (isa(Γ.grid.parent,Gridap.Geometry.GridView))
+      Ltrian=view(Ltrian,ifacet_to_facet)
+    end
+    lhₑ_dofs = get_cell_dof_ids(trial_hybridizable, Ltrian)
     lhₑ_dofs = lazy_map(m, lhₑ_dofs.args[skeleton_fields]...)
 
-    Ω = trial_hybridizable[first(bulk_fields)]
-    Ω = get_triangulation(Ω)
-
     m = Gridap.Fields.BlockMap(length(bulk_fields), bulk_fields)
-    uhph_dofs = get_cell_dof_ids(trial_hybridizable, Ω)
+    uhph_dofs = get_cell_dof_ids(trial_hybridizable, Γ)
   # This last step is needed as get_cell_dof_ids(...) returns as many blocks
   # as fields in trial, regardless of the FEspaces defined on Ω or not
     uhph_dofs = lazy_map(m, uhph_dofs.args[bulk_fields]...)
@@ -188,9 +204,10 @@ function convert_cell_wise_dofs_array_to_facet_dofs_array(
       cells_around_facets,
       cell_wise_facets,
       cell_dofs::AbstractVector{<:AbstractVector{<:Real}},
-      facet_dofs_ids::AbstractVector{<:AbstractVector{<:Integer}})
-    glue = _generate_glue_among_facet_and_cell_wise_dofs_arrays(
-   cells_around_facets, cell_wise_facets, facet_dofs_ids)
+      facet_dofs_ids::AbstractVector{<:AbstractVector{<:Integer}},
+      ifacet_to_facet::AbstractVector{<:Integer}=IdentityVector(length(facet_dofs_ids)))
+   glue = _generate_glue_among_facet_and_cell_wise_dofs_arrays(
+   cells_around_facets, cell_wise_facets, facet_dofs_ids, ifacet_to_facet)
     k = ExtractFacetDofsFromCellDofsMap(cell_dofs)
     [lazy_map(k, glue)]
 end
@@ -198,30 +215,36 @@ end
 function _generate_glue_among_facet_and_cell_wise_dofs_arrays(
   cells_around_facets,
   cell_wise_facets,
-  facet_dof_ids::AbstractVector{<:AbstractVector{<:Integer}})
-
-    c1 = array_cache(cells_around_facets)
-    c2 = array_cache(cell_wise_facets)
-    c3 = array_cache(facet_dof_ids)
-
-    result = Vector{NTuple{3,Int}}(undef, length(facet_dof_ids))
-    current = 1
-    ndofs = 0
-    for facet_gid = 1:length(cells_around_facets)
-        cell_gid = Gridap.Arrays.getindex!(c1, cells_around_facets, facet_gid)[1]
-        current_cell_facets = Gridap.Arrays.getindex!(c2, cell_wise_facets, cell_gid)
-        pos = 1
-        for facet_gid_in_cell in current_cell_facets
-            ndofs = length(Gridap.Arrays.getindex!(c3, facet_dof_ids, facet_gid_in_cell))
-            if (facet_gid == facet_gid_in_cell)
-                break
-            else
-                pos = pos + ndofs
-            end
-        end
-        result[facet_gid] = (cell_gid, pos, ndofs)
-    end
-    result
+  facet_dof_ids::AbstractVector{<:AbstractVector{<:Integer}},
+  ifacet_to_facet::AbstractVector{<:Integer}=IdentityVector(length(facet_dof_ids)))
+
+  @check length(ifacet_to_facet) == length(facet_dof_ids)
+
+  facet_to_ifacet=Dict([facet=>ifacet for (ifacet,facet) in enumerate(ifacet_to_facet)])
+
+  c1 = array_cache(cells_around_facets)
+  c2 = array_cache(cell_wise_facets)
+  c3 = array_cache(facet_dof_ids)
+
+  result = Vector{NTuple{3,Int}}(undef, length(facet_dof_ids))
+  current = 1
+  ndofs = 0
+  for (ifacet,facet_gid) in enumerate(ifacet_to_facet)
+      cell_gid = Gridap.Arrays.getindex!(c1, cells_around_facets, facet_gid)[1]
+      current_cell_facets = Gridap.Arrays.getindex!(c2, cell_wise_facets, cell_gid)
+      pos = 1
+      for facet_gid_in_cell in current_cell_facets
+          ndofs = length(Gridap.Arrays.getindex!(c3, facet_dof_ids, facet_gid_in_cell))
+
+          if (facet_gid == facet_gid_in_cell)
+              break
+          else
+              pos = pos + ndofs
+          end
+      end
+      result[ifacet] = (cell_gid, pos, ndofs)
+  end
+  result
 end
 
 function convert_cell_wise_dofs_array_to_facet_dofs_array(
@@ -362,10 +385,18 @@ function Gridap.FESpaces.get_cell_fe_data(
   tglue::SkeletonGlue) where Dc
     model = tglue.trian.model
     if Dc == num_cell_dims(model)
-        sface_to_data
+      if (isa(tglue.trian.grid.parent,Gridap.Geometry.GridView))
+        cell_to_parent_cell=tglue.trian.grid.parent.cell_to_parent_cell
+        sface_to_data = lazy_map(Reindex(sface_to_data),cell_to_parent_cell)
+      end
+      sface_to_data
     else
         @assert Dc == num_cell_dims(model) - 1
         cell_wise_facets = _get_cell_wise_facets(model)
+        if (isa(tglue.trian.grid.parent,Gridap.Geometry.GridView))
+          cell_to_parent_cell=tglue.trian.grid.parent.cell_to_parent_cell
+          cell_wise_facets = lazy_map(Reindex(cell_wise_facets),cell_to_parent_cell)
+        end
         restrict_facet_dof_ids_to_cell_boundary(cell_wise_facets, sface_to_data)
     end
 end
@@ -441,7 +472,12 @@ end
 function Gridap.FESpaces.get_cell_fe_data(
   fun::typeof(Gridap.FESpaces.get_cell_is_dirichlet),sface_to_data,sglue::FaceToFaceGlue,tglue::SkeletonGlue)
     model = tglue.trian.model
+    grid  = tglue.trian.grid.parent
     cell_wise_facets = _get_cell_wise_facets(model)
+    if (isa(grid,Gridap.Geometry.GridView))
+      cell_to_parent_cell = grid.cell_to_parent_cell
+      cell_wise_facets=lazy_map(Reindex(cell_wise_facets),cell_to_parent_cell)
+    end
     fdofscb = restrict_facet_dof_ids_to_cell_boundary(cell_wise_facets, sface_to_data)
     _generate_cell_is_dirichlet(fdofscb)
 end
diff --git a/src/Skeleton.jl b/src/Skeleton.jl
index 2ff3f60..aec093a 100644
--- a/src/Skeleton.jl
+++ b/src/Skeleton.jl
@@ -15,44 +15,45 @@ struct SkeletonGrid{Dc,Dp,P,A,B,C,D} <: Grid{Dc,Dp}
   cell_lface_nodes::B
   ftype_freffe::C
   cell_lface_ftype::D
-  function SkeletonGrid(grid::Grid)
-    D = num_cell_dims(grid)
-    Dp = num_point_dims(grid)
-    node_coord = get_node_coordinates(grid)
-    cell_ctype = get_cell_type(grid)
-    ctype_reffe       = get_reffes(grid)
-    @notimplementedif length(ctype_reffe) != 1
-    reffe = first(ctype_reffe)
-    freffes = get_reffaces(ReferenceFE{D-1},reffe)
-    @notimplementedif length(freffes) != 1
-    ftype_freffe = [first(freffes),]
-
-    ####
-    ctype_lface_ftype = map(reffe->get_face_type(reffe,D-1),ctype_reffe)
-    cell_lface_ftype  = expand_cell_data(ctype_lface_ftype,cell_ctype)
-    ###
-
-    cell_nodes = get_cell_node_ids(grid)
-    function f(reffe)
-      lface_to_lnodes=get_face_nodes(reffe,D-1)
-      ArrayBlock(lface_to_lnodes,[true for i=1:length(lface_to_lnodes)])
-    end
-    ctype_lface_to_lnodes = map(f,ctype_reffe)
-    cell_lface_to_lnodes  = expand_cell_data(ctype_lface_to_lnodes,cell_ctype)
-    m=lazy_map(Reindex,cell_nodes)
-    fi = testitem(cell_lface_to_lnodes)
-    mi = testitem(m)
-    T = return_type(mi, fi)
-    cell_lface_nodes=LazyArray(T,m,cell_lface_to_lnodes)
-
-    A = typeof(node_coord)
-    B = typeof(cell_lface_nodes)
-    C = typeof(ftype_freffe)
-    E = typeof(cell_lface_ftype)
-    P = typeof(grid)
-    new{D-1,Dp,P,A,B,C,E}(
-      grid,node_coord,cell_lface_nodes,ftype_freffe,cell_lface_ftype)
+end
+
+function SkeletonGrid(grid::Grid)
+  D = num_cell_dims(grid)
+  Dp = num_point_dims(grid)
+  node_coord = get_node_coordinates(grid)
+  cell_ctype = get_cell_type(grid)
+  ctype_reffe       = get_reffes(grid)
+  @notimplementedif length(ctype_reffe) != 1
+  reffe = first(ctype_reffe)
+  freffes = get_reffaces(ReferenceFE{D-1},reffe)
+  @notimplementedif length(freffes) != 1
+  ftype_freffe = [first(freffes),]
+
+  ####
+  ctype_lface_ftype = map(reffe->get_face_type(reffe,D-1),ctype_reffe)
+  cell_lface_ftype  = expand_cell_data(ctype_lface_ftype,cell_ctype)
+  ###
+
+  cell_nodes = get_cell_node_ids(grid)
+  function f(reffe)
+    lface_to_lnodes=get_face_nodes(reffe,D-1)
+    ArrayBlock(lface_to_lnodes,[true for i=1:length(lface_to_lnodes)])
   end
+  ctype_lface_to_lnodes = map(f,ctype_reffe)
+  cell_lface_to_lnodes  = expand_cell_data(ctype_lface_to_lnodes,cell_ctype)
+  m=lazy_map(Reindex,cell_nodes)
+  fi = testitem(cell_lface_to_lnodes)
+  mi = testitem(m)
+  T = return_type(mi, fi)
+  cell_lface_nodes=LazyArray(T,m,cell_lface_to_lnodes)
+
+  A = typeof(node_coord)
+  B = typeof(cell_lface_nodes)
+  C = typeof(ftype_freffe)
+  E = typeof(cell_lface_ftype)
+  P = typeof(grid)
+  SkeletonGrid{D-1,Dp,P,A,B,C,E}(
+    grid,node_coord,cell_lface_nodes,ftype_freffe,cell_lface_ftype)
 end
 
 Geometry.get_node_coordinates(a::SkeletonGrid) = a.node_coord
@@ -98,7 +99,7 @@ function Gridap.Arrays.return_type(k::Reindex,i::VectorBlock{<:Vector{<:Integer}
   TR=eltype(k.values)
   VectorBlock{Vector{TR}}
 end
-# Default return_values fails because one(::VectorBlock{<:Vector{<:Integer}}) is NOT defined
+# Default return_value fails because one(::VectorBlock{<:Vector{<:Integer}}) is NOT defined
 # This justifies why I had to define the following function
 function Gridap.Arrays.return_value(k::Reindex,i::VectorBlock{<:Vector{<:Integer}})
   evaluate(k,i)
@@ -115,32 +116,107 @@ struct SkeletonTriangulation{Dc,Dp,A,B,C} <: Triangulation{Dc,Dp}
   grid::B
   sign_flip::C
   glue::Gridap.Geometry.FaceToCellGlue
-  function SkeletonTriangulation(model::DiscreteModel)
-    A     = typeof(model)
-    D     = num_cell_dims(model)
-    mgrid = get_grid(model)
-    fgrid = Grid(ReferenceFE{D-1},model)
-    glue =  Gridap.Geometry.FaceToCellGlue(get_grid_topology(model),
-                                           mgrid,
-                                           fgrid,
-                                           collect(1:num_facets(model)),
-                                           Fill(Int8(1),num_facets(model)))
-    sgrid = SkeletonGrid(mgrid)
-    B = typeof(sgrid)
-
-    # Generate sign_flip
-    # TO-DO: here I am reusing the machinery for global RT FE spaces.
-    #        Sure there is a way to decouple this from global RT FE spaces.
-    function _get_sign_flip(model)
-      basis,reffe_args,reffe_kwargs = ReferenceFE(raviart_thomas,Float64,0)
-      cell_reffe = ReferenceFE(model,basis,reffe_args...;reffe_kwargs...)
-      Gridap.FESpaces.get_sign_flip(model,cell_reffe)
-    end
-    sign_flip=_get_sign_flip(model)
+end
 
-    C = typeof(sign_flip)
-    new{D-1,D,A,B,C}(model,sgrid,sign_flip,glue)
-  end
+# Generate sign_flip
+# TO-DO: here I am reusing the machinery for global RT FE spaces.
+#        Sure there is a way to decouple this from global RT FE spaces.
+function _get_sign_flip(model)
+  basis,reffe_args,reffe_kwargs = ReferenceFE(raviart_thomas,Float64,0)
+  cell_reffe = ReferenceFE(model,basis,reffe_args...;reffe_kwargs...)
+  Gridap.FESpaces.get_sign_flip(model,cell_reffe)
+end
+
+function SkeletonTriangulation(model::DiscreteModel,
+                               sign_flip=_get_sign_flip(model))
+  A     = typeof(model)
+  D     = num_cell_dims(model)
+  mgrid = get_grid(model)
+  fgrid = Grid(ReferenceFE{D-1},model)
+  glue =  Gridap.Geometry.FaceToCellGlue(get_grid_topology(model),
+                                         mgrid,
+                                         fgrid,
+                                         collect(1:num_facets(model)),
+                                         Fill(Int8(1),num_facets(model)))
+  sgrid = SkeletonGrid(mgrid)
+  B = typeof(sgrid)
+
+  C = typeof(sign_flip)
+  SkeletonTriangulation{D-1,D,A,B,C}(model,sgrid,sign_flip,glue)
+end
+
+function _find_faces_touched_by_cells(cell_to_parent_cell,
+                                      cells_to_facets::Table,
+                                      facets_to_cells::Table)
+
+    nfacets = 0
+    touched=Dict{Int,Bool}()
+    cache=array_cache(cell_to_parent_cell)
+    for i in eachindex(cell_to_parent_cell)
+      cell = getindex!(cache,cell_to_parent_cell,i)
+      a = cells_to_facets.ptrs[cell]
+      b = cells_to_facets.ptrs[cell+1]-1
+      for j=a:b
+         facet=cells_to_facets.data[j]
+         if !(facet in keys(touched))
+            touched[facet]=true
+            nfacets+=1
+         end
+      end
+    end
+    touched=Dict{Int,Bool}()
+    T = eltype(eltype(cells_to_facets))
+    ifacet_to_facet = zeros(T,nfacets)
+    facet_to_icell = ones(T,length(facets_to_cells))
+    nfacets = 0
+    for i in eachindex(cell_to_parent_cell)
+      cell = getindex!(cache,cell_to_parent_cell,i)
+      a = cells_to_facets.ptrs[cell]
+      b = cells_to_facets.ptrs[cell+1]-1
+      for j=a:b
+         facet=cells_to_facets.data[j]
+         if !(facet in keys(touched))
+            touched[facet]=true
+            nfacets+=1
+            ifacet_to_facet[nfacets]=facet
+            c = facets_to_cells.ptrs[facet]
+            d = facets_to_cells.ptrs[facet+1]
+            if ((d-c)==2)
+              if (facets_to_cells.data[d-1]==cell)
+                facet_to_icell[facet]=2
+              end
+            end
+         end
+      end
+    end
+    ifacet_to_facet, facet_to_icell
+end
+
+function SkeletonTriangulation(model::DiscreteModel{D},
+                               cell_to_parent_cell::AbstractVector{<:Integer},
+                               sign_flip=_get_sign_flip(model)) where D
+  topo            = get_grid_topology(model)
+  cells_to_facets = Table(get_faces(topo,D,D-1))
+  facets_to_cells = Table(get_faces(topo,D-1,D))
+  ifacet_to_facet, facet_to_icell = _find_faces_touched_by_cells(cell_to_parent_cell,
+                                                                 cells_to_facets,
+                                                                 facets_to_cells)
+  facet_grid      = Grid(ReferenceFE{D-1},model)
+  cell_grid       = get_grid(model)
+  ifacet_grid     = view(facet_grid,ifacet_to_facet)
+
+  glue =  Gridap.Geometry.FaceToCellGlue(topo,
+                                         cell_grid,
+                                         ifacet_grid,
+                                         ifacet_to_facet,
+                                         facet_to_icell)
+  A = typeof(model)
+  cell_grid = view(cell_grid,cell_to_parent_cell)
+  sgrid = SkeletonGrid(cell_grid)
+  B     = typeof(sgrid)
+  sign_flip=lazy_map(Reindex(sign_flip),cell_to_parent_cell)
+  C = typeof(sign_flip)
+  SkeletonTriangulation{D-1,D,A,B,C}(model,sgrid,sign_flip,glue)
 end
 
 Skeleton(args...) = SkeletonTriangulation(args...)
@@ -201,8 +277,28 @@ function Gridap.FESpaces.get_cell_fe_data(
   sface_to_data
 end
 
-function Geometry.best_target(a::BodyFittedTriangulation{Dca},
-                              b::BodyFittedTriangulation{Dcb}) where {Dca,Dcb}
+function Geometry.best_target(a::BodyFittedTriangulation{1},
+                              b::BodyFittedTriangulation{2}) where {Dca,Dcb}
+  _best_target(a,b)
+end
+
+function Geometry.best_target(a::BodyFittedTriangulation{2},
+                              b::BodyFittedTriangulation{1}) where {Dca,Dcb}
+  _best_target(a,b)
+end
+
+function Geometry.best_target(a::BodyFittedTriangulation{2},
+                              b::BodyFittedTriangulation{3}) where {Dca,Dcb}
+  _best_target(a,b)
+end
+
+function Geometry.best_target(a::BodyFittedTriangulation{3},
+                              b::BodyFittedTriangulation{2}) where {Dca,Dcb}
+  _best_target(a,b)
+end
+
+function _best_target(a::BodyFittedTriangulation{Dca},
+                      b::BodyFittedTriangulation{Dcb}) where {Dca,Dcb}
   @assert Dca==Dcb-1 || Dca-1==Dcb
   @assert get_background_model(a)===get_background_model(b)
   Skeleton(get_background_model(a))
@@ -210,15 +306,56 @@ end
 
 # TO-DO: dirty. I cannot check whether a===b, as a and b might be created from scratch
 #               along the process
+# TO-DO: fix current implementation of this function. It is nothing but a quick
+#        and dirty workaround. Possibly using TriangulationView and GridView.
 function Geometry.best_target(a::SkeletonTriangulation{Dc},
                               b::SkeletonTriangulation{Dc}) where {Dc}
-  a
+  agrid=a.grid.parent
+  bgrid=b.grid.parent
+  if (typeof(agrid)==typeof(bgrid))
+    a
+  elseif (isa(agrid,Gridap.Geometry.GridView))
+    a
+  elseif (isa(bgrid,Gridap.Geometry.GridView))
+    b
+  else
+    @notimplemented
+  end
+
 end
 
+# TO-DO: fix current implementation of this function. It is nothing but a quick
+#        and dirty workaround. Possibly using TriangulationView and GridView.
 function CellData.change_domain_ref_ref(a::CellField,
                                         ttrian::SkeletonTriangulation,
-                                        sglue::SkeletonGlue,tglue::SkeletonGlue)
-  a
+                                        sglue::SkeletonGlue,
+                                        tglue::SkeletonGlue)
+
+  sgrid=sglue.trian.grid.parent
+  tgrid=tglue.trian.grid.parent
+  @assert get_background_model(sglue.trian)===get_background_model(tglue.trian)
+  if (typeof(sgrid)==typeof(tgrid) && !isa(sgrid,Gridap.Geometry.GridView))
+     a
+  else
+     if (!isa(sgrid,Gridap.Geometry.GridView) &&
+          isa(tgrid,Gridap.Geometry.GridView))
+        cell_to_parent_cell=tgrid.cell_to_parent_cell
+        a_field = Gridap.CellData.get_data(a)
+        a_field_view = lazy_map(Reindex(a_field),cell_to_parent_cell)
+        Gridap.CellData.similar_cell_field(a,a_field_view,ttrian,ReferenceDomain())
+     elseif (isa(sgrid,Gridap.Geometry.GridView) &&
+              !isa(tgrid,Gridap.Geometry.GridView))
+        cell_to_parent_cell=sgrid.cell_to_parent_cell
+        a_field = Gridap.CellData.get_data(a)
+        a_field_view = lazy_map(Reindex(a_field),cell_to_parent_cell)
+        Gridap.CellData.similar_cell_field(a,a_field_view,ttrian,ReferenceDomain())
+     else
+        scell_to_parent_cell=sgrid.cell_to_parent_cell
+        tcell_to_parent_cell=tgrid.cell_to_parent_cell
+        @notimplementedif !all(scell_to_parent_cell==tcell_to_parent_cell)
+        a
+     end
+  end
 end
 
 function Geometry.get_glue(trian::SkeletonTriangulation{D},::Val{D}) where D
@@ -249,12 +386,17 @@ function Geometry.get_glue(trian::SkeletonTriangulation{d},::Val{D}) where {d,D}
     poly = get_polytope(reffe)
     num_faces(poly,d)
   end
-  # Avoid allocations here
-  tcell_lface_mface = lazy_map(cell_ctype,cell_cell) do ctype, cell
+  function f(ctype, cell)
     nlfaces = ctype_nlfaces[ctype]
     fill(cell,nlfaces)
   end
-  tcell_lface_mface_map = _setup_tcell_lface_mface_map(d,trian.model,trian.glue)
+  # Avoid allocations here
+  # TO-THINK. I am not 100% sure if tcell_lface_mface now
+  #           holds the data it was planned to hold.
+  #           It does not cause trouble as it turns out tcell_lface_mface
+  #           is not currently being used.
+  tcell_lface_mface = lazy_map(f,cell_ctype,cell_cell)
+  tcell_lface_mface_map = _setup_tcell_lface_mface_map(d,trian.model,trian.grid.parent,trian.glue)
   SkeletonGlue(trian,tcell_lface_mface,tcell_lface_mface_map)
 end
 
@@ -268,38 +410,42 @@ function CellData.change_domain_ref_ref(
 
   if isa(strian,Triangulation{D,D})
     b=_restrict_to_skeleton_cell_field(ttrian.model,
+                                       ttrian.grid.parent,
                                        ttrian.glue,
                                        tglue.tcell_lface_mface_map,
                                        a)
   elseif isa(strian,Triangulation{D-1,D})
-    b=_restrict_to_skeleton_facet_field(ttrian.model,ttrian.glue,a)
+    b=_restrict_to_skeleton_facet_field(ttrian.model,ttrian.grid.parent,ttrian.glue,a)
   end
   CellData.similar_cell_field(a,b,ttrian,ReferenceDomain())
 end
 
 function _restrict_to_skeleton_cell_field(model,
+                                          cell_grid,
                                           glue,
                                           tface_to_mface_map,
                                           cell_fe_basis::Gridap.CellData.CellField)
   D = num_cell_dims(model)
   Gridap.Helpers.@check isa(get_triangulation(cell_fe_basis),Triangulation{D,D})
-  cell_a_q = _transform_cell_to_cell_lface_array(glue,
+  cell_a_q = _transform_cell_to_cell_lface_array(glue,cell_grid,
          Gridap.CellData.get_data(cell_fe_basis);
          add_naive_innermost_block_level=true)
   lazy_map(Broadcasting(∘),cell_a_q,tface_to_mface_map)
 end
 
 function _transform_cell_to_cell_lface_array(glue,
-                                            cell_array::Fill;
-                                            add_naive_innermost_block_level=false)
-    d = Gridap.Arrays.CompressedArray([cell_array.value,],Fill(1,length(cell_array)))
-    _transform_cell_to_cell_lface_array(glue,d;add_naive_innermost_block_level=add_naive_innermost_block_level)
+                                             cell_grid,
+                                             cell_array::Fill;
+                                             add_naive_innermost_block_level=false)
+    d = Gridap.Arrays.CompressedArray([cell_array.value,],Fill(1,num_cells(cell_grid)))
+    _transform_cell_to_cell_lface_array(glue,cell_grid,d;add_naive_innermost_block_level=add_naive_innermost_block_level)
 end
 
 
 function _transform_cell_to_cell_lface_array(glue,
-                                            cell_array::Gridap.Arrays.CompressedArray;
-                                            add_naive_innermost_block_level=false)
+                                             cell_grid,
+                                             cell_array::Gridap.Arrays.CompressedArray;
+                                             add_naive_innermost_block_level=false)
   T=typeof(cell_array.values[1])
   ctype_to_vector_block=
     Vector{Gridap.Fields.VectorBlock{T}}(undef,length(glue.ctype_to_lface_to_ftype))
@@ -316,13 +462,14 @@ function _transform_cell_to_cell_lface_array(glue,
   if add_naive_innermost_block_level
     ctype_to_vector_block=collect(lazy_map(AddNaiveInnerMostBlockLevelMap(),ctype_to_vector_block))
   end
-  Gridap.Arrays.CompressedArray(ctype_to_vector_block,glue.cell_to_ctype)
+  Gridap.Arrays.CompressedArray(ctype_to_vector_block,get_cell_type(cell_grid))
 end
 
 function _transform_cell_to_cell_lface_array(glue,
+  cell_grid,
   cell_array::AbstractVector;
   add_naive_innermost_block_level=false)
-  ctype_to_vector_block=SkeletonVectorFromCellVector(glue,cell_array)
+  ctype_to_vector_block=SkeletonVectorFromCellVector(glue,cell_grid,cell_array)
   if add_naive_innermost_block_level
     ctype_to_vector_block=collect(lazy_map(AddNaiveInnerMostBlockLevelMap(),ctype_to_vector_block))
   end
@@ -339,16 +486,18 @@ function _transform_cell_to_cell_lface_array(
 end
 
 function _restrict_to_skeleton_facet_field(model,
+                                           cell_grid,
                                            glue,
                                            facet_fe_function::Gridap.FESpaces.SingleFieldFEFunction)
   D = num_cell_dims(model)
   Gridap.Helpers.@check isa(get_triangulation(facet_fe_function),Triangulation{D-1,D})
   facet_field_array=Gridap.CellData.get_data(facet_fe_function)
   cell_wise_facets_ids=_get_cell_wise_facets(model)
-  SkeletonVectorFromFacetVector(glue,cell_wise_facets_ids,facet_field_array)
+  SkeletonVectorFromFacetVector(glue,cell_grid,cell_wise_facets_ids,facet_field_array)
 end
 
 function _restrict_to_skeleton_facet_field(model,
+                                           cell_grid,
                                            glue,
                                            facet_fe_basis::Gridap.CellData.CellField)
 
@@ -357,11 +506,13 @@ function _restrict_to_skeleton_facet_field(model,
 
   _transform_face_to_cell_lface_expanded_array(
     glue,
+    cell_grid,
     Gridap.CellData.get_data(facet_fe_basis))
 end
 
 function _transform_face_to_cell_lface_expanded_array(
   glue,
+  cell_grid,
   face_array::Gridap.Arrays.LazyArray{<:Fill{typeof(transpose)}})
   Gridap.Helpers.@check typeof(face_array.args[1]) <: Gridap.Arrays.CompressedArray
 
@@ -371,12 +522,13 @@ function _transform_face_to_cell_lface_expanded_array(
     v[i]=evaluate(transpose,face_array.args[1].values[i])
   end
   face_array_compressed=Gridap.Arrays.CompressedArray(v,face_array.args[1].ptrs)
-  _transform_face_to_cell_lface_expanded_array(glue,face_array_compressed)
+  _transform_face_to_cell_lface_expanded_array(glue,cell_grid,face_array_compressed)
 end
 
 
 function _transform_face_to_cell_lface_expanded_array(glue,
-                                                     face_array::Fill)
+                                                      cell_grid,
+                                                      face_array::Fill)
   T=typeof(face_array.value)
   ctype_to_vector_block=
     Vector{Gridap.Fields.VectorBlock{T}}(undef,length(glue.ctype_to_lface_to_ftype))
@@ -391,11 +543,12 @@ function _transform_face_to_cell_lface_expanded_array(glue,
      end
      ctype_to_vector_block[ctype]=Gridap.Fields.ArrayBlock(v,t)
   end
-  Gridap.Arrays.CompressedArray(ctype_to_vector_block,glue.cell_to_ctype)
+  Gridap.Arrays.CompressedArray(ctype_to_vector_block,get_cell_type(cell_grid))
 end
 
 function _transform_face_to_cell_lface_expanded_array(glue,
-                                                     face_array::Gridap.Arrays.CompressedArray)
+                                                      cell_grid,
+                                                      face_array::Gridap.Arrays.CompressedArray)
   ftype_to_block_layout=_get_block_layout(face_array.values)
   T=eltype(face_array.values[1])
   if length(ftype_to_block_layout[1][1]) == 1
@@ -449,7 +602,7 @@ function _transform_face_to_cell_lface_expanded_array(glue,
      end
      ctype_to_vector_block[ctype]=Gridap.Fields.ArrayBlock(vf1,tf1)
   end
-  Gridap.Arrays.CompressedArray(ctype_to_vector_block,glue.cell_to_ctype)
+  Gridap.Arrays.CompressedArray(ctype_to_vector_block,get_cell_type(cell_grid))
 end
 
 function _get_block_layout(fields_array::AbstractArray{<:AbstractArray{<:Gridap.Fields.Field}})
@@ -463,7 +616,7 @@ end
 
 function CellData.change_domain_phys_phys(
   a::CellField,ttrian::SkeletonTriangulation,sglue::FaceToFaceGlue,tglue::SkeletonGlue)
-  sface_to_field = get_data(a)
+  sface_to_field = Gridap.CellData.get_data(a)
   mface_to_sface = sglue.mface_to_tface
   tcell_lface_mface = tglue.tcell_lface_mface
   mface_to_field = extend(sface_to_field,mface_to_sface)
@@ -479,13 +632,12 @@ Returns a cell-wise array which, for each cell, and each facet within the cell,
 returns the unit outward normal to the boundary of the cell.
 """
 function get_cell_normal_vector(s::SkeletonTriangulation)
-  cell_lface_normal=_get_cell_normal_vector(s.model, s.glue, _cell_lface_to_nref)
+  cell_lface_normal=_get_cell_normal_vector(s.model, s.grid.parent, s.glue, _cell_lface_to_nref)
   GenericCellField(cell_lface_normal,s,ReferenceDomain())
 end
 
 function _cell_lface_to_nref(args...)
-  model,glue = args[1],first(args[2:end])
-  cell_grid = get_grid(model)
+  model, cell_grid, glue = args[1],args[2],first(args[3:end])
   ## Reference normal
   function f(r)
     p = Gridap.ReferenceFEs.get_polytope(r)
@@ -496,7 +648,7 @@ function _cell_lface_to_nref(args...)
     lface_pindex_to_n
   end
   ctype_lface_pindex_to_nref = map(f, get_reffes(cell_grid))
-  SkeletonCompressedVector(ctype_lface_pindex_to_nref,glue)
+  SkeletonCompressedVector(ctype_lface_pindex_to_nref,cell_grid,glue)
 end
 
 """
@@ -506,6 +658,7 @@ returns the unit outward normal to the boundary of the cell owner of the facet.
 function get_cell_owner_normal_vector(s::SkeletonTriangulation)
   cell_owner_lface_normal=_get_cell_normal_vector(
      s.model,
+     s.grid.parent,
      s.glue,
      _cell_lface_to_owner_nref,
      s.sign_flip)
@@ -513,25 +666,28 @@ function get_cell_owner_normal_vector(s::SkeletonTriangulation)
 end
 
 function _cell_lface_to_owner_nref(args...)
-  model,glue,sign_flip = args
-  cell_lface_to_nref=_cell_lface_to_nref(model,glue)
+  model,cell_grid,glue,sign_flip = args
+  cell_lface_to_nref=_cell_lface_to_nref(model,cell_grid,glue)
   SkeletonOwnerNref(cell_lface_to_nref,sign_flip)
 end
 
-function _get_cell_normal_vector(model,glue,cell_lface_to_nref::Function,sign_flip=nothing)
-  cell_grid = get_grid(model)
+function _get_cell_normal_vector(model,
+                                 cell_grid,
+                                 glue,
+                                 cell_lface_to_nref::Function,
+                                 sign_flip=nothing)
 
-  cell_lface_to_nref = cell_lface_to_nref(model,glue,sign_flip)
+  cell_lface_to_nref = cell_lface_to_nref(model,cell_grid,glue,sign_flip)
   cell_lface_s_nref = lazy_map(Gridap.Fields.constant_field,cell_lface_to_nref)
 
   # Inverse of the Jacobian transpose
   cell_q_x = get_cell_map(cell_grid)
   cell_q_Jt = lazy_map(∇,cell_q_x)
   cell_q_invJt = lazy_map(Operation(Gridap.Fields.pinvJt),cell_q_Jt)
-  cell_lface_q_invJt = _transform_cell_to_cell_lface_array(glue, cell_q_invJt)
+  cell_lface_q_invJt = _transform_cell_to_cell_lface_array(glue, cell_grid, cell_q_invJt)
 
   # Change of domain
-  cell_lface_s_q = _setup_tcell_lface_mface_map(num_cell_dims(model)-1,model,glue)
+  cell_lface_s_q = _setup_tcell_lface_mface_map(num_cell_dims(model)-1,model,cell_grid,glue)
 
   cell_lface_s_invJt = lazy_map(∘,cell_lface_q_invJt,cell_lface_s_q)
   #face_s_n =
@@ -541,16 +697,17 @@ function _get_cell_normal_vector(model,glue,cell_lface_to_nref::Function,sign_fl
   #Fields.MemoArray(face_s_n)
 end
 
-function _setup_tcell_lface_mface_map(d,model,glue)
-  ctype_to_lface_to_pindex_to_qcoords=Gridap.Geometry._compute_face_to_q_vertex_coords_body(d,model,glue)
-  cell_lface_to_q_vertex_coords = SkeletonCompressedVector(
-                                  ctype_to_lface_to_pindex_to_qcoords.ctype_lface_pindex_to_value,
-                                  glue)
+function _setup_tcell_lface_mface_map(d,model,cell_grid,glue)
+  ctype_to_lface_to_pindex_to_qcoords=
+     Gridap.Geometry._compute_face_to_q_vertex_coords_body(d,model,glue)
+  cell_lface_to_q_vertex_coords =
+    SkeletonCompressedVector(ctype_to_lface_to_pindex_to_qcoords.ctype_lface_pindex_to_value,
+                             cell_grid,
+                             glue)
   f(p) = Gridap.ReferenceFEs.get_shapefuns(
          Gridap.ReferenceFEs.LagrangianRefFE(Float64,Gridap.ReferenceFEs.get_polytope(p),1))
 
   ################ TO-IMPROVE
-  cell_grid = get_grid(model)
   D = num_cell_dims(model)
   ctype_reffe = get_reffes(cell_grid)
   Gridap.Helpers.@notimplementedif length(ctype_reffe) != 1
@@ -562,13 +719,14 @@ function _setup_tcell_lface_mface_map(d,model,glue)
 
   ftype_to_shapefuns = map( f,  ftrian_reffes)
   face_to_shapefuns = expand_cell_data(ftype_to_shapefuns,glue.face_to_ftype)
-  cell_to_lface_to_shapefuns = transform_face_to_cell_lface_array(glue,face_to_shapefuns)
+  cell_to_lface_to_shapefuns = transform_face_to_cell_lface_array(glue,cell_grid,face_to_shapefuns)
   lazy_map(Gridap.Fields.linear_combination,
            cell_lface_to_q_vertex_coords,
            cell_to_lface_to_shapefuns)
 end
 
 function transform_face_to_cell_lface_array(glue,
+                                            cell_grid,
                                             face_array::Gridap.Arrays.CompressedArray,
                                             f::Function=identity)
   T=typeof(f(face_array.values[1]))
@@ -585,7 +743,7 @@ function transform_face_to_cell_lface_array(glue,
      end
      ctype_to_vector_block[ctype]=Gridap.Fields.ArrayBlock(v,t)
   end
-  Gridap.Arrays.CompressedArray(ctype_to_vector_block,glue.cell_to_ctype)
+  Gridap.Arrays.CompressedArray(ctype_to_vector_block,get_cell_type(cell_grid))
 end
 
 function Gridap.ReferenceFEs.expand_cell_data(
diff --git a/src/SkeletonArrays.jl b/src/SkeletonArrays.jl
index bafbcf1..bc5bf8d 100644
--- a/src/SkeletonArrays.jl
+++ b/src/SkeletonArrays.jl
@@ -1,9 +1,12 @@
-struct SkeletonCompressedVector{T,G<:Gridap.Geometry.FaceToCellGlue} <: AbstractVector{Gridap.Fields.VectorBlock{T}}
+struct SkeletonCompressedVector{T,
+                                R,
+                                G<:Gridap.Geometry.FaceToCellGlue} <: AbstractVector{Gridap.Fields.VectorBlock{T}}
   ctype_lface_pindex_to_value::Vector{Vector{Vector{T}}}
+  cell_grid::R
   glue::G
 end
 
-function _compressed_vector_from_glue(::Type{T}, glue) where T
+function _compressed_vector_from_glue(::Type{T}, glue, cell_grid) where T
   ctype_to_vector_block=
     Vector{Gridap.Fields.VectorBlock{T}}(undef,length(glue.ctype_to_lface_to_ftype))
   for ctype=1:length(glue.ctype_to_lface_to_ftype)
@@ -13,11 +16,11 @@ function _compressed_vector_from_glue(::Type{T}, glue) where T
      t.=true
      ctype_to_vector_block[ctype]=Gridap.Fields.ArrayBlock(v,t)
   end
-  Gridap.Arrays.CompressedArray(ctype_to_vector_block, glue.cell_to_ctype)
+  Gridap.Arrays.CompressedArray(ctype_to_vector_block, get_cell_type(cell_grid))
 end
 
 function Gridap.Arrays.array_cache(a::SkeletonCompressedVector{T}) where {T}
-   _compressed_vector_from_glue(T,a.glue)
+   _compressed_vector_from_glue(T,a.glue,a.cell_grid)
 end
 
 function Base.getindex(a::SkeletonCompressedVector,cell::Integer)
@@ -27,7 +30,7 @@ end
 
 function Gridap.Arrays.getindex!(cache,a::SkeletonCompressedVector,cell::Integer)
    vb=cache[cell]
-   ctype=a.glue.cell_to_ctype[cell]
+   ctype=get_cell_type(a.cell_grid)[cell]
    for lface=1:length(vb)
     p = a.glue.cell_to_lface_to_pindex.ptrs[cell]-1
     pindex = a.glue.cell_to_lface_to_pindex.data[p+lface]
@@ -36,13 +39,13 @@ function Gridap.Arrays.getindex!(cache,a::SkeletonCompressedVector,cell::Integer
    vb
 end
 
-Base.size(a::SkeletonCompressedVector) = (length(a.glue.cell_to_ctype),)
+Base.size(a::SkeletonCompressedVector) = (num_cells(a.cell_grid),)
 
 Base.IndexStyle(::Type{<:SkeletonCompressedVector}) = IndexLinear()
 
 function Gridap.Arrays.lazy_map(k::Gridap.Fields.LinearCombinationMap,
                   ::Type{T},b::SkeletonCompressedVector,c::Fill) where T
-  d = Gridap.Arrays.CompressedArray([c.value,],Fill(1,length(c)))
+  d = Gridap.Arrays.CompressedArray([c.value,],Fill(1,num_cells(b.cell_grid)))
   lazy_map(k,T,b,d)
 end
 
@@ -50,7 +53,7 @@ function Gridap.Arrays.lazy_map(k::Gridap.Fields.LinearCombinationMap,
                                 ::Type{T},
                                 b::SkeletonCompressedVector,
                                 c::Gridap.Arrays.CompressedArray{<:Gridap.Fields.VectorBlock}) where T
-  if c.ptrs === b.glue.cell_to_ctype || c.ptrs == b.glue.cell_to_ctype
+  if c.ptrs === get_cell_type(b.cell_grid) || c.ptrs == get_cell_type(b.cell_grid)
     ET=eltype(T)
     ctype_lface_pindex_to_r = Vector{Vector{Vector{ET}}}(undef,length(b.ctype_lface_pindex_to_value))
     for (ctype, lface_pindex_to_value) in enumerate(b.ctype_lface_pindex_to_value)
@@ -67,16 +70,12 @@ function Gridap.Arrays.lazy_map(k::Gridap.Fields.LinearCombinationMap,
       end
       ctype_lface_pindex_to_r[ctype] = lface_pindex_to_r
     end
-    SkeletonCompressedVector(ctype_lface_pindex_to_r,b.glue)
+    SkeletonCompressedVector(ctype_lface_pindex_to_r,b.cell_grid,b.glue)
   else
     Gridap.Helpers.@notimplemented
   end
 end
 
-
-
-
-
 function Gridap.Arrays.lazy_map(
   k::typeof(Gridap.Arrays.evaluate),
   ::Type{T},
@@ -99,7 +98,7 @@ function Gridap.Arrays.lazy_map(
     end
     ctype_lface_pindex_to_r[ctype] = lface_pindex_to_r
   end
-  SkeletonCompressedVector(ctype_lface_pindex_to_r,b.glue)
+  SkeletonCompressedVector(ctype_lface_pindex_to_r,b.cell_grid,b.glue)
 end
 
 function Gridap.Arrays.lazy_map(
@@ -125,7 +124,7 @@ function Gridap.Arrays.lazy_map(
     end
     ctype_lface_pindex_to_r[ctype] = lface_pindex_to_r
   end
-  SkeletonCompressedVector(ctype_lface_pindex_to_r,b.glue)
+  SkeletonCompressedVector(ctype_lface_pindex_to_r,b.cell_grid,b.glue)
 end
 
 function Gridap.Arrays.lazy_map(
@@ -150,33 +149,35 @@ function Gridap.Arrays.lazy_map(
     end
     ctype_lface_pindex_to_r[ctype] = lface_pindex_to_r
   end
-  SkeletonCompressedVector(ctype_lface_pindex_to_r,b.glue)
+  SkeletonCompressedVector(ctype_lface_pindex_to_r,b.cell_grid,b.glue)
 end
 
 
 
-struct SkeletonVectorFromFacetVector{T,G,C,V} <: AbstractVector{Gridap.Fields.VectorBlock{T}}
+struct SkeletonVectorFromFacetVector{T,G,F,C,V} <: AbstractVector{Gridap.Fields.VectorBlock{T}}
   glue::G
+  cell_grid::F
   cell_wise_facets_ids::C
   facet_vector::V
-  function SkeletonVectorFromFacetVector(glue,cell_wise_facets_ids,facet_vector)
+  function SkeletonVectorFromFacetVector(glue,cell_grid,cell_wise_facets_ids,facet_vector)
     G=typeof(glue)
+    F=typeof(cell_grid)
     C=typeof(cell_wise_facets_ids)
     V=typeof(facet_vector)
     T=eltype(V)
-    new{T,G,C,V}(glue,cell_wise_facets_ids,facet_vector)
+    new{T,G,F,C,V}(glue,cell_grid,cell_wise_facets_ids,facet_vector)
   end
 end
 
-Base.size(a::SkeletonVectorFromFacetVector) = (length(a.glue.cell_to_ctype),)
+Base.size(a::SkeletonVectorFromFacetVector) = (length(get_cell_type(a.cell_grid)),)
 
 Base.IndexStyle(::Type{<:SkeletonVectorFromFacetVector}) = IndexLinear()
 
 function Gridap.Arrays.array_cache(a::SkeletonVectorFromFacetVector{T}) where T
    fvc=array_cache(a.facet_vector)
    cwfc=array_cache(a.cell_wise_facets_ids)
-   vbc=_compressed_vector_from_glue(T,a.glue)
-   fvc=_compressed_vector_from_glue(typeof(fvc),a.glue)
+   vbc=_compressed_vector_from_glue(T,a.glue,a.cell_grid)
+   fvc=_compressed_vector_from_glue(typeof(fvc),a.glue,a.cell_grid)
    for ctype=1:length(a.glue.ctype_to_lface_to_ftype)
     num_facets=length(a.glue.ctype_to_lface_to_ftype[ctype])
     for lface=1:num_facets
@@ -240,7 +241,7 @@ function Gridap.Arrays.lazy_map(
   af=_cell_lfacet_vector_to_facet_vector(b.glue,a)
   bf=b.facet_vector
   bfx=lazy_map(evaluate,bf,af)
-  SkeletonVectorFromFacetVector(b.glue,b.cell_wise_facets_ids,bfx)
+  SkeletonVectorFromFacetVector(b.glue,b.cell_grid,b.cell_wise_facets_ids,bfx)
 end
 
 function _cell_lfacet_vector_to_facet_vector(
@@ -276,28 +277,30 @@ function Gridap.Arrays.lazy_map(
   ::typeof(∇),
   b::SkeletonVectorFromFacetVector)
   ∇bf=lazy_map(∇,b.cell_vector)
-  SkeletonVectorFromFacetVector(b.glue,b.cell_wise_facets_ids,∇bf)
+  SkeletonVectorFromFacetVector(b.glue,b.cell_grid,b.cell_wise_facets_ids,∇bf)
 end
 
 
-struct SkeletonVectorFromCellVector{T,G,V} <: AbstractVector{Gridap.Fields.VectorBlock{T}}
+struct SkeletonVectorFromCellVector{T,G,R,V} <: AbstractVector{Gridap.Fields.VectorBlock{T}}
   glue::G
+  cell_grid::R
   cell_vector::V
-  function SkeletonVectorFromCellVector(glue,cell_vector)
+  function SkeletonVectorFromCellVector(glue,cell_grid,cell_vector)
     G=typeof(glue)
+    R=typeof(cell_grid)
     V=typeof(cell_vector)
     T=eltype(V)
-    new{T,G,V}(glue,cell_vector)
+    new{T,G,R,V}(glue,cell_grid,cell_vector)
   end
 end
 
-Base.size(a::SkeletonVectorFromCellVector) = (length(a.glue.cell_to_ctype),)
+Base.size(a::SkeletonVectorFromCellVector) = (num_cells(a.cell_grid),)
 
 Base.IndexStyle(::Type{<:SkeletonVectorFromCellVector}) = IndexLinear()
 
 function Gridap.Arrays.array_cache(a::SkeletonVectorFromCellVector{T}) where T
    cvc=array_cache(a.cell_vector)
-   vbc=_compressed_vector_from_glue(T,a.glue)
+   vbc=_compressed_vector_from_glue(T,a.glue,a.cell_grid)
    cvc,vbc
 end
 
diff --git a/test/Distributed/DarcyHDGTests.jl b/test/Distributed/DarcyHDGTests.jl
new file mode 100644
index 0000000..8581eba
--- /dev/null
+++ b/test/Distributed/DarcyHDGTests.jl
@@ -0,0 +1,138 @@
+module DarcyHDGTests
+
+using Test
+using Gridap
+using GridapDistributed
+using PartitionedArrays
+using Gridap
+using FillArrays
+using Gridap.Geometry
+using GridapHybrid
+
+function main(parts)
+  partition = (0,1,0,1,0,1)
+  cells = (2,2,2)
+  model = CartesianDiscreteModel(parts,partition,cells)
+  order=1
+  solve_darcy_lhdg(model,order)
+end
+
+u2(x) = VectorValue(1+x[1],1+x[2])
+Gridap.divergence(::typeof(u2)) = (x) -> 2
+p(x) = -3.14
+∇p2(x) = VectorValue(0,0,0)
+Gridap.∇(::typeof(p)) = ∇p
+f2(x) = u2(x) + ∇p2(x)
+# Normal component of u(x) on Neumann boundary
+function g2(x)
+  tol=1.0e-14
+  if (abs(x[2])<tol)
+    return -x[2] #-x[1]-x[2]
+  elseif (abs(x[2]-1.0)<tol)
+    return x[2] # x[1]+x[2]
+  end
+  Gridap.Helpers.@check false
+end
+
+
+#3D problem
+u3(x) = VectorValue(1+x[1],1+x[2],1+x[3])
+Gridap.divergence(::typeof(u3)) = (x) -> 3
+∇p3(x) = VectorValue(0,0,0)
+f3(x) = u3(x) + ∇p3(x)
+function g3(x) # Normal component of u(x) on Neumann boundary
+  @assert false
+end
+
+function ufg(D::Int)
+   if (D==2)
+    u2,f2,g2
+   elseif (D==3)
+    u3,f3,g3
+   end
+end
+
+function dirichlet_tags(D::Int)
+  if (D==2)
+    collect(5:8)
+  elseif (D==3)
+    collect(21:26)
+  end
+end
+
+
+
+function solve_darcy_lhdg(model,order)
+    # Geometry
+    D = num_cell_dims(model)
+
+    dtags=dirichlet_tags(D)
+    u,f,_ = ufg(D)
+
+    Ω = Triangulation(ReferenceFE{D},model)
+    Γ = Triangulation(ReferenceFE{D-1},model)
+    ∂K = GridapHybrid.Skeleton(model)
+
+    # Reference FEs
+    order  = 1
+    reffeᵤ = ReferenceFE(lagrangian,VectorValue{D,Float64},order;space=:P)
+    reffeₚ = ReferenceFE(lagrangian,Float64,order-1;space=:P)
+    reffeₗ = ReferenceFE(lagrangian,Float64,order;space=:P)
+
+    # Define test FESpaces
+    V = TestFESpace(Ω  , reffeᵤ; conformity=:L2)
+    Q = TestFESpace(Ω  , reffeₚ; conformity=:L2)
+    M = TestFESpace(Γ,
+                    reffeₗ;
+                    conformity=:L2,
+                    dirichlet_tags=dtags)
+    Y = MultiFieldFESpace([V,Q,M])
+
+    # Define trial FEspaces
+    U = TrialFESpace(V)
+    P = TrialFESpace(Q)
+    L = TrialFESpace(M,p)
+    X = MultiFieldFESpace([U, P, L])
+
+    # FE formulation params
+    τ = 1.0 # HDG stab parameter
+
+    degree = 2*order+1
+    dΩ     = Measure(Ω,degree)
+    n      = get_cell_normal_vector(∂K)
+    nₒ     = get_cell_owner_normal_vector(∂K)
+    d∂K    = Measure(∂K,degree)
+
+    # (uh,ph,lh) = xh
+    # (vh,qh,mh) = yh
+
+    # ∫( vh⋅uh - (∇⋅vh)*ph - ∇(qh)⋅uh )dΩ
+    # (vh⋅n)
+    # (vh⋅n)*lh
+    # ∫((vh⋅n)*lh)d∂K
+    # ∫(qh*(uh⋅n))d∂K
+    # ∫(τ*qh*ph*(n⋅nₒ))d∂K
+    # ∫(mh*(uh⋅n))d∂K
+    # ∫(τ*mh*ph*(n⋅nₒ))d∂K
+    # ∫(τ*mh*lh*(n⋅nₒ))d∂K
+
+    a((uh,ph,lh),(vh,qh,mh)) = ∫( vh⋅uh - (∇⋅vh)*ph - ∇(qh)⋅uh )dΩ +
+                              ∫((vh⋅n)*lh)d∂K +
+                              #∫(qh*(uh⋅n+τ*(ph-lh)*n⋅no))*d∂K
+                              ∫(qh*(uh⋅n))d∂K +
+                              ∫(τ*qh*ph*(n⋅nₒ))d∂K -
+                              ∫(τ*qh*lh*(n⋅nₒ))d∂K +
+                              #∫(mh*(uh⋅n+τ*(ph-lh)*n⋅no))*d∂K
+                              ∫(mh*(uh⋅n))d∂K +
+                              ∫(τ*mh*ph*(n⋅nₒ))d∂K -
+                              ∫(τ*mh*lh*(n⋅nₒ))d∂K
+    l((vh,qh,mh)) = ∫( vh⋅f + qh*(∇⋅u))*dΩ
+
+    op=HybridAffineFEOperator((u,v)->(a(u,v),l(v)), X, Y, [1,2], [3])
+    xh=solve(op)
+
+    uh,_=xh
+    e = u -uh
+    @test sqrt(sum(∫(e⋅e)dΩ)) < 1.0e-12
+  end
+end # module
diff --git a/test/Distributed/mpi/runtests.jl b/test/Distributed/mpi/runtests.jl
new file mode 100644
index 0000000..0ad0018
--- /dev/null
+++ b/test/Distributed/mpi/runtests.jl
@@ -0,0 +1,31 @@
+module MPITests
+
+using MPI
+using Test
+
+#Sysimage
+sysimage=nothing
+if length(ARGS)==1
+   @assert isfile(ARGS[1]) "$(ARGS[1]) must be a valid Julia sysimage file"
+   sysimage=ARGS[1]
+end
+
+mpidir = @__DIR__
+testdir = joinpath(mpidir,"..")
+repodir = joinpath(mpidir,"../../..")
+function run_driver(procs,file,sysimage)
+  mpiexec() do cmd
+    if sysimage!=nothing
+       extra_args="-J$(sysimage)"
+       run(`$cmd -n $procs $(Base.julia_cmd()) $(extra_args) --project=$repodir $(joinpath(mpidir,file))`)
+    else
+       run(`$cmd -n $procs $(Base.julia_cmd()) --project=$repodir $(joinpath(mpidir,file))`)
+    end
+    @test true
+  end
+end
+
+run_driver(1,"runtests_body.jl",sysimage) # Check that the degenerated case works
+run_driver(4,"runtests_body.jl",sysimage)
+
+end # module
diff --git a/test/Distributed/mpi/runtests_body.jl b/test/Distributed/mpi/runtests_body.jl
new file mode 100644
index 0000000..1213f71
--- /dev/null
+++ b/test/Distributed/mpi/runtests_body.jl
@@ -0,0 +1,35 @@
+module NP4
+# All test running on 4 procs here
+
+using PartitionedArrays
+const PArrays = PartitionedArrays
+using MPI
+using GridapPETSc
+
+include("../DarcyHDGTests.jl")
+
+if ! MPI.Initialized()
+  MPI.Init()
+end
+
+function all_tests(parts)
+  display(parts)
+  t = PArrays.PTimer(parts,verbose=true)
+  PArrays.tic!(t)
+  options = "-ksp_type cg -pc_type gamg -ksp_monitor"
+  GridapPETSc.with(args=split(options)) do
+     DarcyHDGTests.main(parts)
+  end
+  PArrays.toc!(t,"DarcyHDGTests")
+  display(t)
+end
+
+if MPI.Comm_size(MPI.COMM_WORLD) == 8
+  prun(all_tests,mpi,(2,2,2))
+elseif MPI.Comm_size(MPI.COMM_WORLD) == 1
+  prun(all_tests,mpi,(1,1,1))
+else
+  MPI.Abort(MPI.COMM_WORLD,0)
+end
+
+end #module
diff --git a/test/Distributed/sequential/DarcyHDGTests.jl b/test/Distributed/sequential/DarcyHDGTests.jl
new file mode 100644
index 0000000..d58bf73
--- /dev/null
+++ b/test/Distributed/sequential/DarcyHDGTests.jl
@@ -0,0 +1,6 @@
+module DarcyHDGTestsSeq
+using PartitionedArrays
+include("../DarcyHDGTests.jl")
+prun(DarcyHDGTests.main,sequential,(1,1,1))
+prun(DarcyHDGTests.main,sequential,(2,2,2))
+end # module
diff --git a/test/runtests.jl b/test/runtests.jl
index aa32229..0650572 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,4 +9,5 @@ module Tests
   @testset "DarcyHDGTests" begin include("DarcyHDGTests.jl") end
   @testset "LinearElasticityHDGTests" begin include("LinearElasticityHDGTests.jl") end
   @testset "MultiFieldLagrangeMultipliersTests" begin include("MultiFieldLagrangeMultipliersTests.jl") end
+  @testset "Distributed/DarcyHDGTests" begin include("Distributed/sequential/DarcyHDGTests.jl") end
 end # module