diff --git a/Manifest.toml b/Manifest.toml
index 4849931..38789ed 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.9.2"
 manifest_format = "2.0"
-project_hash = "3cb3efc1daa833e06773d949d12011bd0cc455e9"
+project_hash = "547bf31d5dc2e6452ad979d2c526bfa03b483ccd"
 
 [[deps.AbstractFFTs]]
 deps = ["LinearAlgebra"]
@@ -25,9 +25,9 @@ version = "0.4.4"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "76289dc51920fdc6e0013c872ba9551d54961c24"
+git-tree-sha1 = "f8c724a2066b2d37d0234fe4022ec67987022d00"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "3.6.2"
+version = "4.0.0"
 weakdeps = ["StaticArrays"]
 
     [deps.Adapt.extensions]
@@ -50,9 +50,9 @@ version = "1.1.1"
 
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra", "Requires", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "f83ec24f76d4c8f525099b2ac475fc098138ec31"
+git-tree-sha1 = "bbec08a37f8722786d87bedf84eae19c020c4efa"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.4.11"
+version = "7.7.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
@@ -72,9 +72,9 @@ version = "7.4.11"
 
 [[deps.ArrayLayouts]]
 deps = ["FillArrays", "LinearAlgebra"]
-git-tree-sha1 = "9a731850434825d183af39c6e6cd0a1c32dd7e20"
+git-tree-sha1 = "b08a4043e1c14096ef8efe4dd97e07de5cacf240"
 uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
-version = "1.4.2"
+version = "1.4.5"
 weakdeps = ["SparseArrays"]
 
     [deps.ArrayLayouts.extensions]
@@ -85,9 +85,9 @@ uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 
 [[deps.AutoHashEquals]]
 deps = ["Pkg"]
-git-tree-sha1 = "16d2cc66d4de2b4cac1b0942671e07edbfb922c5"
+git-tree-sha1 = "daaeb6f7f77b88c072a83a2451801818acb5c63b"
 uuid = "15f4f7f2-30c1-5605-9d31-71845cf9641f"
-version = "2.0.0"
+version = "2.1.0"
 
 [[deps.BSON]]
 git-tree-sha1 = "2208958832d6e1b59e49f53697483a84ca8d664e"
@@ -99,9 +99,9 @@ uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
 [[deps.BlockArrays]]
 deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra"]
-git-tree-sha1 = "54cd829dd26330c42e1cf9df68470dd4df602c61"
+git-tree-sha1 = "fc69cbdb4277042f72c6e59cbc7024fbe3034b89"
 uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-version = "0.16.38"
+version = "0.16.39"
 
 [[deps.Bzip2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -122,9 +122,9 @@ version = "1.16.1+1"
 
 [[deps.CircularArrays]]
 deps = ["OffsetArrays"]
-git-tree-sha1 = "61bc114e595167090b4cbcb7305ddeacd4274f16"
+git-tree-sha1 = "3f7b8a37359ae592cfa7aca7f811da045deff222"
 uuid = "7a955b69-7140-5f4e-a0ed-f168c5e2e749"
-version = "1.3.2"
+version = "1.3.3"
 
 [[deps.CodecZlib]]
 deps = ["TranscodingStreams", "Zlib_jll"]
@@ -157,9 +157,9 @@ version = "0.3.0"
 
 [[deps.Compat]]
 deps = ["UUIDs"]
-git-tree-sha1 = "8a62af3e248a8c4bad6b32cbbe663ae02275e32c"
+git-tree-sha1 = "886826d76ea9e72b35fcd000e535588f7b60f21d"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.10.0"
+version = "4.10.1"
 weakdeps = ["Dates", "LinearAlgebra"]
 
     [deps.Compat.extensions]
@@ -208,9 +208,9 @@ version = "1.15.1"
 
 [[deps.Distances]]
 deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
-git-tree-sha1 = "5225c965635d8c21168e32a12954675e7bea1151"
+git-tree-sha1 = "66c4c81f259586e8f002eacebc177e1fb06363b0"
 uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
-version = "0.10.10"
+version = "0.10.11"
 
     [deps.Distances.extensions]
     DistancesChainRulesCoreExt = "ChainRulesCore"
@@ -249,9 +249,9 @@ version = "2.5.0+0"
 
 [[deps.FFTW]]
 deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
-git-tree-sha1 = "b4fbdd20c889804969571cc589900803edda16b7"
+git-tree-sha1 = "ec22cbbcd01cba8f41eecd7d44aac1f23ee985e3"
 uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-version = "1.7.1"
+version = "1.7.2"
 
 [[deps.FFTW_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -282,15 +282,20 @@ uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
 [[deps.FillArrays]]
 deps = ["LinearAlgebra", "Random"]
-git-tree-sha1 = "a20eaa3ad64254c61eeb5f230d9306e937405434"
+git-tree-sha1 = "5b93957f6dcd33fc343044af3d48c215be2562f1"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "1.6.1"
-weakdeps = ["SparseArrays", "Statistics"]
+version = "1.9.3"
 
     [deps.FillArrays.extensions]
+    FillArraysPDMatsExt = "PDMats"
     FillArraysSparseArraysExt = "SparseArrays"
     FillArraysStatisticsExt = "Statistics"
 
+    [deps.FillArrays.weakdeps]
+    PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+    Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
 [[deps.FiniteDiff]]
 deps = ["ArrayInterface", "LinearAlgebra", "Requires", "Setfield", "SparseArrays"]
 git-tree-sha1 = "c6e4a1fbe73b31a3dea94b1da449503b8830c306"
@@ -375,19 +380,19 @@ version = "1.3.1"
 
 [[deps.Gridap]]
 deps = ["AbstractTrees", "BSON", "BlockArrays", "Combinatorics", "DataStructures", "DocStringExtensions", "FastGaussQuadrature", "FileIO", "FillArrays", "ForwardDiff", "JLD2", "JSON", "LineSearches", "LinearAlgebra", "NLsolve", "NearestNeighbors", "PolynomialBases", "QuadGK", "Random", "SparseArrays", "SparseMatricesCSR", "StaticArrays", "Test", "WriteVTK"]
-git-tree-sha1 = "666d4d66201e7cc9142fe940f8e42ac6fb616291"
+git-tree-sha1 = "0ff991a6beee6cbe9e0804ef3cb368cbe4ed83ab"
 repo-rev = "master"
 repo-url = "https://github.com/gridap/Gridap.jl.git"
 uuid = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
-version = "0.17.20"
+version = "0.17.21"
 
 [[deps.GridapDistributed]]
 deps = ["BlockArrays", "FillArrays", "Gridap", "LinearAlgebra", "MPI", "PartitionedArrays", "SparseArrays", "SparseMatricesCSR", "WriteVTK"]
-git-tree-sha1 = "a61a30f052a9a9d2a1548c80d1a130e78c62cb3c"
+git-tree-sha1 = "acad4231797aa33b8ad3f636ad8e8d075eee024a"
 repo-rev = "master"
 repo-url = "https://github.com/gridap/GridapDistributed.jl.git"
 uuid = "f9701e48-63b3-45aa-9a63-9bc6c271f355"
-version = "0.3.1"
+version = "0.3.5"
 
 [[deps.GridapGmsh]]
 deps = ["Gridap", "GridapDistributed", "Libdl", "Metis", "PartitionedArrays", "gmsh_jll"]
@@ -397,22 +402,24 @@ version = "0.7.0"
 
 [[deps.GridapP4est]]
 deps = ["ArgParse", "FillArrays", "Gridap", "GridapDistributed", "Libdl", "MPI", "P4est_wrapper", "PartitionedArrays", "Test"]
-git-tree-sha1 = "8188694b39e7afafd19539d9864db96ad67159b9"
-repo-rev = "uniform-3d"
+git-tree-sha1 = "db6555cabb76b1362f0044f5a7b90ef145c86389"
+repo-rev = "main"
 repo-url = "https://github.com/gridap/GridapP4est.jl.git"
 uuid = "c2c8e14b-f5fd-423d-9666-1dd9ad120af9"
-version = "0.3.1"
+version = "0.3.4"
 
 [[deps.GridapPETSc]]
 deps = ["Gridap", "GridapDistributed", "Libdl", "LinearAlgebra", "MPI", "PETSc_jll", "PartitionedArrays", "Random", "SparseArrays", "SparseMatricesCSR"]
-git-tree-sha1 = "fc15fcc0833736ce264b0454edad40e8f3176b7f"
+git-tree-sha1 = "fdef4c7834299334044da32dae4003d3020ebf5f"
+repo-rev = "add_missing_api"
+repo-url = "https://github.com/gridap/GridapPETSc.jl.git"
 uuid = "bcdc36c2-0c3e-11ea-095a-c9dadae499f1"
 version = "0.5.0"
 
 [[deps.GridapSolvers]]
-deps = ["AbstractTrees", "ArgParse", "BlockArrays", "FillArrays", "Gridap", "GridapDistributed", "GridapP4est", "GridapPETSc", "IterativeSolvers", "LinearAlgebra", "MPI", "PartitionedArrays", "Printf", "SparseArrays", "SparseMatricesCSR"]
-git-tree-sha1 = "f2176c38ef5bf05e5e10603770c36416b01a2bae"
-repo-rev = "develop"
+deps = ["AbstractTrees", "ArgParse", "BlockArrays", "FillArrays", "Gridap", "GridapDistributed", "GridapP4est", "GridapPETSc", "IterativeSolvers", "LinearAlgebra", "MPI", "P4est_wrapper", "PartitionedArrays", "Printf", "SparseArrays", "SparseMatricesCSR"]
+git-tree-sha1 = "e4808b9bfa48625c23c9f7a2494a108fb464e99c"
+repo-rev = "block-solvers"
 repo-url = "https://github.com/gridap/GridapSolvers.jl.git"
 uuid = "6d3209ee-5e3c-4db7-a716-942eb12ed534"
 version = "0.2.0"
@@ -429,10 +436,10 @@ uuid = "9b13fd28-a010-5f03-acff-a1bbcff69959"
 version = "1.0.0"
 
 [[deps.IntelOpenMP_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "ad37c091f7d7daf900963171600d7c1c5c3ede32"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "5fdf2fe6724d8caabf43b557b84ce53f3b7e2f6b"
 uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
-version = "2023.2.0+0"
+version = "2024.0.2+0"
 
 [[deps.InteractiveUtils]]
 deps = ["Markdown"]
@@ -445,15 +452,15 @@ version = "0.2.2"
 
 [[deps.IterativeSolvers]]
 deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays"]
-git-tree-sha1 = "1169632f425f79429f245113b775a0e3d121457c"
+git-tree-sha1 = "b435d190ef8369cf4d79cc9dd5fba88ba0165307"
 uuid = "42fd0dbc-a981-5370-80f2-aaf504508153"
-version = "0.9.2"
+version = "0.9.3"
 
 [[deps.JLD2]]
-deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "Pkg", "Printf", "Reexport", "Requires", "TranscodingStreams", "UUIDs"]
-git-tree-sha1 = "c11d691a0dc8e90acfa4740d293ade57f68bfdbb"
+deps = ["FileIO", "MacroTools", "Mmap", "OrderedCollections", "Pkg", "PrecompileTools", "Printf", "Reexport", "Requires", "TranscodingStreams", "UUIDs"]
+git-tree-sha1 = "c2d0f45afcb5f6209155670bffd100c3b4937ea3"
 uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-version = "0.4.35"
+version = "0.4.40"
 
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
@@ -469,15 +476,15 @@ version = "0.21.4"
 
 [[deps.JpegTurbo_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "6f2675ef130a300a112286de91973805fcc5ffbc"
+git-tree-sha1 = "60b1194df0a3298f460063de985eae7b01bc011a"
 uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
-version = "2.1.91+0"
+version = "3.0.1+0"
 
 [[deps.LLVMOpenMP_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "f689897ccbe049adb19a065c495e75f372ecd42b"
+deps = ["Artifacts", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "d986ce2d884d49126836ea94ed5bfb0f12679713"
 uuid = "1d63c593-3942-5779-bab2-d838dc0a180e"
-version = "15.0.4+0"
+version = "15.0.7+0"
 
 [[deps.LZO_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -561,9 +568,9 @@ version = "2.36.0+0"
 
 [[deps.LightXML]]
 deps = ["Libdl", "XML2_jll"]
-git-tree-sha1 = "e129d9391168c677cd4800f5c0abb1ed8cb3794f"
+git-tree-sha1 = "3a994404d3f6709610701c7dabfc03fed87a81f8"
 uuid = "9c8b4983-aa76-5018-a973-4c85ecc9e179"
-version = "0.9.0"
+version = "0.9.1"
 
 [[deps.LineSearches]]
 deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"]
@@ -607,10 +614,10 @@ uuid = "d00139f3-1899-568f-a2f0-47f597d42d70"
 version = "5.1.2+0"
 
 [[deps.MKL_jll]]
-deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
-git-tree-sha1 = "eb006abbd7041c28e0d16260e50a24f8f9104913"
+deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl"]
+git-tree-sha1 = "72dc3cf284559eb8f53aa593fe62cb33f83ed0c0"
 uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
-version = "2023.2.0+0"
+version = "2024.0.0+0"
 
 [[deps.MMG_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "LinearElasticity_jll", "Pkg", "SCOTCH_jll"]
@@ -620,9 +627,9 @@ version = "5.6.0+0"
 
 [[deps.MPI]]
 deps = ["Distributed", "DocStringExtensions", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "PkgVersion", "PrecompileTools", "Requires", "Serialization", "Sockets"]
-git-tree-sha1 = "b4d8707e42b693720b54f0b3434abee6dd4d947a"
+git-tree-sha1 = "4e3136db3735924f96632a5b40a5979f1f53fa07"
 uuid = "da04e1cc-30fd-572f-bb4f-1f8673147195"
-version = "0.20.16"
+version = "0.20.19"
 
     [deps.MPI.extensions]
     AMDGPUExt = "AMDGPU"
@@ -634,27 +641,27 @@ version = "0.20.16"
 
 [[deps.MPICH_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
-git-tree-sha1 = "8a5b4d2220377d1ece13f49438d71ad20cf1ba83"
+git-tree-sha1 = "2ee75365ca243c1a39d467e35ffd3d4d32eef11e"
 uuid = "7cb0a576-ebde-5e09-9194-50597f1243b4"
-version = "4.1.2+0"
+version = "4.1.2+1"
 
 [[deps.MPIPreferences]]
 deps = ["Libdl", "Preferences"]
-git-tree-sha1 = "781916a2ebf2841467cda03b6f1af43e23839d85"
+git-tree-sha1 = "8f6af051b9e8ec597fa09d8885ed79fd582f33c9"
 uuid = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267"
-version = "0.1.9"
+version = "0.1.10"
 
 [[deps.MPItrampoline_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML"]
-git-tree-sha1 = "6979eccb6a9edbbb62681e158443e79ecc0d056a"
+git-tree-sha1 = "8eeb3c73bbc0ca203d0dc8dad4008350bbe5797b"
 uuid = "f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748"
-version = "5.3.1+0"
+version = "5.3.1+1"
 
 [[deps.MacroTools]]
 deps = ["Markdown", "Random"]
-git-tree-sha1 = "9ee1618cbf5240e6d4e0371d6f24065083f60c48"
+git-tree-sha1 = "b211c553c199c111d998ecdaf7623d1b89b69f93"
 uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
-version = "0.5.11"
+version = "0.5.12"
 
 [[deps.Markdown]]
 deps = ["Base64"]
@@ -683,9 +690,9 @@ version = "1.4.0"
 
 [[deps.MicrosoftMPI_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "a7023883872e52bc29bcaac74f19adf39347d2d5"
+git-tree-sha1 = "b01beb91d20b0d1312a9471a36017b5b339d26de"
 uuid = "9237b28f-5490-5468-be7b-bb81f5f5e6cf"
-version = "10.1.4+0"
+version = "10.1.4+1"
 
 [[deps.Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
@@ -714,9 +721,9 @@ version = "1.0.2"
 
 [[deps.NearestNeighbors]]
 deps = ["Distances", "StaticArrays"]
-git-tree-sha1 = "2c3726ceb3388917602169bed973dbc97f1b51a8"
+git-tree-sha1 = "ded64ff6d4fdd1cb68dfcbb818c69e144a5b2e4c"
 uuid = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
-version = "0.4.13"
+version = "0.4.16"
 
 [[deps.NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
@@ -729,10 +736,13 @@ uuid = "baad4e97-8daa-5946-aac2-2edac59d34e1"
 version = "7.6.2+2"
 
 [[deps.OffsetArrays]]
-deps = ["Adapt"]
-git-tree-sha1 = "2ac17d29c523ce1cd38e27785a7d23024853a4bb"
+git-tree-sha1 = "6a731f2b5c03157418a20c12195eb4b74c8f8621"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.12.10"
+version = "1.13.0"
+weakdeps = ["Adapt"]
+
+    [deps.OffsetArrays.extensions]
+    OffsetArraysAdaptExt = "Adapt"
 
 [[deps.OpenBLAS32_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
@@ -769,9 +779,9 @@ uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
 version = "0.5.5+0"
 
 [[deps.OrderedCollections]]
-git-tree-sha1 = "2e73fe17cac3c62ad1aebe70d44c963c3cfdc3e3"
+git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.6.2"
+version = "1.6.3"
 
 [[deps.P4est_jll]]
 deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "MPICH_jll", "MPIPreferences", "MPItrampoline_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "Pkg", "TOML", "Zlib_jll"]
@@ -781,7 +791,9 @@ version = "2.8.1+2"
 
 [[deps.P4est_wrapper]]
 deps = ["CEnum", "Libdl", "MPI", "P4est_jll"]
-git-tree-sha1 = "c182e067e6bac213c6e75e4139d01ebfb7ea0a58"
+git-tree-sha1 = "d1dc7da59aefe29ff874970dc141fab07378e791"
+repo-rev = "p4est-2.3+"
+repo-url = "https://github.com/gridap/P4est_wrapper.jl.git"
 uuid = "3743d7c0-8adf-11ea-380b-7d33b0ecc1da"
 version = "0.2.0"
 
@@ -797,10 +809,10 @@ uuid = "8fa3689e-f0b9-5420-9873-adf6ccf46f2d"
 version = "3.15.2+0"
 
 [[deps.PackageCompiler]]
-deps = ["Artifacts", "Glob", "LazyArtifacts", "Libdl", "Pkg", "Printf", "RelocatableFolders", "TOML", "UUIDs"]
-git-tree-sha1 = "d24f0920777c64c0265e9725eb4a3ec450bc36dd"
+deps = ["Artifacts", "Glob", "LazyArtifacts", "Libdl", "Pkg", "Printf", "RelocatableFolders", "TOML", "UUIDs", "p7zip_jll"]
+git-tree-sha1 = "8b880733c61c4a99a069302390de2d057ee166ce"
 uuid = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
-version = "2.1.10"
+version = "2.1.16"
 
 [[deps.Parameters]]
 deps = ["OrderedCollections", "UnPack"]
@@ -810,9 +822,9 @@ version = "0.12.3"
 
 [[deps.Parsers]]
 deps = ["Dates", "PrecompileTools", "UUIDs"]
-git-tree-sha1 = "716e24b21538abc91f6205fd1d8363f39b442851"
+git-tree-sha1 = "8489905bcdbcfac64d1daa51ca07c0d8f0283821"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.7.2"
+version = "2.8.1"
 
 [[deps.PartitionedArrays]]
 deps = ["CircularArrays", "Distances", "FillArrays", "IterativeSolvers", "LinearAlgebra", "MPI", "Printf", "Random", "SparseArrays", "SparseMatricesCSR"]
@@ -912,9 +924,9 @@ version = "0.7.0"
 
 [[deps.Scratch]]
 deps = ["Dates"]
-git-tree-sha1 = "30449ee12237627992a99d5e30ae63e4d78cd24a"
+git-tree-sha1 = "3bac05bc7e74a75fd9cba4295cde4045d9fe2386"
 uuid = "6c6a2e73-6563-6170-7368-637461726353"
-version = "1.2.0"
+version = "1.2.1"
 
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -956,15 +968,19 @@ version = "2.3.1"
     ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 
 [[deps.StaticArrays]]
-deps = ["LinearAlgebra", "Random", "StaticArraysCore"]
-git-tree-sha1 = "0adf069a2a490c47273727e029371b31d44b72b2"
+deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
+git-tree-sha1 = "fba11dbe2562eecdfcac49a05246af09ee64d055"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.6.5"
-weakdeps = ["Statistics"]
+version = "1.8.1"
 
     [deps.StaticArrays.extensions]
+    StaticArraysChainRulesCoreExt = "ChainRulesCore"
     StaticArraysStatisticsExt = "Statistics"
 
+    [deps.StaticArrays.weakdeps]
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
 [[deps.StaticArraysCore]]
 git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d"
 uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
@@ -1034,15 +1050,15 @@ version = "1.0.1"
 
 [[deps.WriteVTK]]
 deps = ["Base64", "CodecZlib", "FillArrays", "LightXML", "TranscodingStreams", "VTKBase"]
-git-tree-sha1 = "7b46936613e41cfe1c6a5897d243ddcab8feabec"
+git-tree-sha1 = "41f0dc2a8f6fd860c266b91fd5cdf4fead65ae69"
 uuid = "64499a7a-5c06-52f2-abe2-ccb03c286192"
-version = "1.18.0"
+version = "1.18.1"
 
 [[deps.XML2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"]
-git-tree-sha1 = "24b81b59bd35b3c42ab84fa589086e19be919916"
+git-tree-sha1 = "801cbe47eae69adc50f36c3caec4758d2650741b"
 uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
-version = "2.11.5+0"
+version = "2.12.2+0"
 
 [[deps.XSLT_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
@@ -1133,10 +1149,10 @@ uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
 version = "5.8.0+0"
 
 [[deps.libpng_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"]
+git-tree-sha1 = "93284c28274d9e75218a416c65ec49d0e0fcdf3d"
 uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
-version = "1.6.38+0"
+version = "1.6.40+0"
 
 [[deps.nghttp2_jll]]
 deps = ["Artifacts", "Libdl"]
diff --git a/Project.toml b/Project.toml
index 8a47ff6..1db192c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,6 +17,7 @@ GridapPETSc = "bcdc36c2-0c3e-11ea-095a-c9dadae499f1"
 GridapSolvers = "6d3209ee-5e3c-4db7-a716-942eb12ed534"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
+P4est_wrapper = "3743d7c0-8adf-11ea-380b-7d33b0ecc1da"
 PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
 PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
 Profile = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
diff --git a/src/FEOperators.jl b/src/FEOperators.jl
deleted file mode 100644
index 085bf2e..0000000
--- a/src/FEOperators.jl
+++ /dev/null
@@ -1,52 +0,0 @@
-"""
-
-"""
-struct FEOperatorMHD <: Gridap.FESpaces.FEOperator
-  op_global :: FEOperator # global operator
-  op_u      :: FEOperator # operator for u
-end
-
-Gridap.FESpaces.get_test(op::FEOperatorMHD)  = get_test(op.op_global)
-Gridap.FESpaces.get_trial(op::FEOperatorMHD) = get_trial(op.op_global)
-
-# Delegate most of the methods to the global operator
-Gridap.FESpaces.allocate_residual(op::FEOperatorMHD,uh) = allocate_residual(op.op_global,uh)
-Gridap.FESpaces.residual!(b::AbstractVector,op::FEOperatorMHD,uh) = residual!(b,op.op_global,uh)
-Gridap.FESpaces.residual_and_jacobian(op::FEOperatorMHD,uh) = residual_and_jacobian(op.op_global,uh)
-Gridap.FESpaces.residual_and_jacobian!(b::AbstractVector,A::AbstractMatrix,op::FEOperatorMHD,uh) = residual_and_jacobian!(b,A,op.op_global,uh)
-
-# When updating the matrix, only update the u-u block
-Gridap.FESpaces.allocate_jacobian(op::FEOperatorMHD,uh) = allocate_jacobian(op.op_global,uh)
-Gridap.FESpaces.jacobian(op::FEOperatorMHD,uh) = jacobian(op.op_global,uh)
-function Gridap.FESpaces.jacobian!(A::AbstractBlockMatrix,op::FEOperatorMHD,xh)
-  uh, ph, jh, φh = xh
-  jacobian!(A[Block(1,1)],op.op_u,uh)
-  return A
-end
-
-function weakform_uu(params,k)
-  fluid = params[:fluid]
-  Ωf  = _interior(params[:model],fluid[:domain])
-  dΩf = Measure(Ωf,2*k)
-
-  α  = fluid[:α]
-  β  = fluid[:β]
-  f  = fluid[:f]
-  ζ  = params[:ζ]
-
-  #function a(u,dv)
-  #  r = ∫(β*(∇(u)⊙∇(dv))) * dΩf
-  #  if ζ !== nothing 
-  #    r += ∫(ζ*(∇⋅u)*(∇⋅dv)) * dΩf
-  #  end
-  #  return r
-  #end
-  a(u,dv) = ∫(β*(∇(u)⊙∇(dv))) * dΩf
-  l(dv) = ∫( dv⋅f ) * dΩf
-  c(u,dv) = ∫( α*dv⋅(conv∘(u,∇(u))) ) * dΩf
-  dc(u,du,dv) = ∫( α*dv⋅( (conv∘(u,∇(du))) + (conv∘(du,∇(u))) ) ) * dΩf
-
-  res(u,dv) = c(u,dv) + a(u,dv) - l(dv)
-  jac(u,du,dv) = dc(u,du,dv) + a(du,dv)
-  return res, jac
-end
diff --git a/src/GridapMHD.jl b/src/GridapMHD.jl
index cd98305..d54e551 100644
--- a/src/GridapMHD.jl
+++ b/src/GridapMHD.jl
@@ -31,7 +31,7 @@ using GridapPETSc.PETSC
 
 using GridapSolvers
 using GridapSolvers.LinearSolvers
-using GridapSolvers.LinearSolvers: allocate_col_vector, allocate_row_vector
+using GridapSolvers.LinearSolvers: allocate_in_domain, allocate_in_range
 
 using PartitionedArrays
 using PartitionedArrays: getany
@@ -45,7 +45,6 @@ include("Li2019/Li2019.jl")
 # Main driver
 include("Fixes.jl")
 include("ExtraFunctions.jl")
-include("FEOperators.jl")
 include("parameters.jl")
 include("Main.jl")
 
diff --git a/src/Li2019/Li2019.jl b/src/Li2019/Li2019.jl
index 18fd7a6..43775ab 100644
--- a/src/Li2019/Li2019.jl
+++ b/src/Li2019/Li2019.jl
@@ -25,10 +25,9 @@ using GridapPETSc
 
 using GridapSolvers
 using GridapSolvers.LinearSolvers
-using GridapSolvers.LinearSolvers: allocate_col_vector, allocate_row_vector
+using GridapSolvers.LinearSolvers: allocate_in_domain, allocate_in_range
 
 include("NLSolvers.jl")
-include("PETScLinearSolvers.jl")
 
 include("blocks.jl")
 include("preconditioner.jl")
diff --git a/src/Li2019/NLSolvers.jl b/src/Li2019/NLSolvers.jl
index 75c6788..4c8e5c1 100644
--- a/src/Li2019/NLSolvers.jl
+++ b/src/Li2019/NLSolvers.jl
@@ -2,7 +2,7 @@
 function Gridap.Algebra.solve!(x::AbstractVector,nls::NewtonRaphsonSolver,op::NonlinearOperator,cache::Nothing)
   b  = residual(op, x)
   A  = jacobian(op, x)
-  dx = allocate_col_vector(A)
+  dx = allocate_in_domain(A)
   ns = numerical_setup(symbolic_setup(nls.ls, A), A)
 
   Gridap.Algebra._solve_nr!(x,A,b,dx,ns,nls,op)
diff --git a/src/Li2019/PETScLinearSolvers.jl b/src/Li2019/PETScLinearSolvers.jl
deleted file mode 100644
index 8302afd..0000000
--- a/src/Li2019/PETScLinearSolvers.jl
+++ /dev/null
@@ -1,67 +0,0 @@
-
-using GridapPETSc: PETScMatrix, PETScVector, PETScLinearSolverSS
-using GridapPETSc.PETSC
-using MPI
-
-mutable struct CachedPETScLinearSolverNS <: NumericalSetup
-  B::PETScMatrix
-  ksp::Ref{KSP}
-  initialized::Bool
-  caches
-  function CachedPETScLinearSolverNS(B::PETScMatrix,caches)
-    new(B,Ref{KSP}(),false,caches)
-  end
-end
-
-function Init(a::CachedPETScLinearSolverNS)
-  @assert Threads.threadid() == 1
-  GridapPETSc._NREFS[] += 1
-  a.initialized = true
-  finalizer(Finalize,a)
-end
-
-function Finalize(ns::CachedPETScLinearSolverNS)
-  if ns.initialized && GridapPETSc.Initialized()
-    if ns.B.comm == MPI.COMM_SELF
-      @check_error_code PETSC.KSPDestroy(ns.ksp)
-    else
-      @check_error_code PETSC.PetscObjectRegisterDestroy(ns.ksp[].ptr)
-    end
-    ns.initialized = false
-    @assert Threads.threadid() == 1
-    GridapPETSc._NREFS[] -= 1
-  end
-  nothing
-end
-
-function get_petsc_caches(A::AbstractMatrix)
-  X = convert(PETScVector,allocate_col_vector(A))
-  B = convert(PETScVector,allocate_col_vector(A))
-  return X,B
-end
-
-function Algebra.numerical_setup(ss::PETScLinearSolverSS,A::AbstractMatrix)
-  B = convert(PETScMatrix,A)
-  caches = get_petsc_caches(A)
-  ns = CachedPETScLinearSolverNS(B,caches)
-  @check_error_code PETSC.KSPCreate(B.comm,ns.ksp)
-  @check_error_code PETSC.KSPSetOperators(ns.ksp[],ns.B.mat[],ns.B.mat[])
-  ss.solver.setup(ns.ksp)
-  @check_error_code PETSC.KSPSetUp(ns.ksp[])
-  Init(ns)
-end
-
-function Algebra.solve!(x::AbstractVector{PetscScalar},ns::CachedPETScLinearSolverNS,b::AbstractVector{PetscScalar})
-  X, B = ns.caches
-  copy!(B,b)
-  @check_error_code PETSC.KSPSolve(ns.ksp[],B.vec[],X.vec[])
-  copy!(x,X)
-  return x
-end
-
-function Algebra.numerical_setup!(ns::CachedPETScLinearSolverNS,A::AbstractMatrix)
-  copy!(ns.B,A)
-  @check_error_code PETSC.KSPSetOperators(ns.ksp[],ns.B.mat[],ns.B.mat[])
-  @check_error_code PETSC.KSPSetUp(ns.ksp[])
-  return ns
-end
diff --git a/src/Li2019/preconditioner.jl b/src/Li2019/preconditioner.jl
index 474b177..4ef416f 100644
--- a/src/Li2019/preconditioner.jl
+++ b/src/Li2019/preconditioner.jl
@@ -25,10 +25,10 @@ mutable struct LI2019_NS <: Gridap.Algebra.NumericalSetup
 end
 
 function allocate_caches(solver::Li2019_Preconditioner,A::AbstractBlockMatrix)
-  du = allocate_col_vector(A[Block(1,1)])
-  dp = allocate_col_vector(A[Block(2,2)])
-  dj = allocate_col_vector(A[Block(3,3)])
-  dφ = allocate_col_vector(A[Block(4,4)])
+  du = allocate_in_domain(A[Block(1,1)])
+  dp = allocate_in_domain(A[Block(2,2)])
+  dj = allocate_in_domain(A[Block(3,3)])
+  dφ = allocate_in_domain(A[Block(4,4)])
   return du, dp, dj, dφ
 end
 
diff --git a/src/Li2019/solver.jl b/src/Li2019/solver.jl
index 7383aad..4824243 100644
--- a/src/Li2019/solver.jl
+++ b/src/Li2019/solver.jl
@@ -155,7 +155,7 @@ function test_preconditioner(op,P)
   # Test preconditioner
   i_am_main(parts) && println(repeat("=",50))
   ns = numerical_setup(symbolic_setup(P,A),A)
-  x = allocate_col_vector(A)
+  x = allocate_in_domain(A)
   solve!(x,ns,b)
   e = norm(x)
   i_am_main(parts) && println("P error: ",e)
@@ -168,8 +168,8 @@ function test_block_solvers(parts,block_solvers,block_mats)
 
   function test_solver(s,m)
     ns = numerical_setup(symbolic_setup(s,m),m)
-    x = allocate_col_vector(m)
-    b = allocate_col_vector(m)
+    x = allocate_in_domain(m)
+    b = allocate_in_domain(m)
     fill!(b,1.0)
     x = solve!(x,ns,b)
     return norm(b - m*x)
diff --git a/src/Main.jl b/src/Main.jl
index 145676e..8f9c5d7 100644
--- a/src/Main.jl
+++ b/src/Main.jl
@@ -291,8 +291,13 @@ end
 function weak_form(params,k)
 
   fluid = params[:fluid]
-  Ωf  = _interior(params[:model],fluid[:domain])
-  dΩf = Measure(Ωf,2*k)
+  Ωf    = _interior(params[:model],fluid[:domain])
+  dΩf   = Measure(Ωf,2*k)
+
+  α, β, γ, σf = fluid[:α], fluid[:β], fluid[:γ], fluid[:σ]
+  f = fluid[:f]
+  B = fluid[:B]
+  ζ = params[:ζ]
 
   solid = params[:solid]
   if solid !== nothing
@@ -301,14 +306,6 @@ function weak_form(params,k)
     σs  = solid[:σ]
   end
 
-  α  = fluid[:α]
-  β  = fluid[:β]
-  γ  = fluid[:γ]
-  f  = fluid[:f]
-  B  = fluid[:B]
-  σf = fluid[:σ]
-  ζ  = params[:ζ]
-
   bcs = params[:bcs]
 
   params_φ = []
@@ -466,6 +463,3 @@ function a_B(x,dy,γ,B,dΩ)
   v_u, v_p, v_j, v_φ = dy
   ∫( -(γ*(j×B)⋅v_u) - (u×B)⋅v_j )*dΩ
 end
-
-
-
diff --git a/test/dev/Dj_solver.jl b/test/dev/Dj_solver.jl
index 825bd15..e61ddf1 100644
--- a/test/dev/Dj_solver.jl
+++ b/test/dev/Dj_solver.jl
@@ -9,9 +9,9 @@ using Gridap.ReferenceFEs
 function test_smoother(s,D_j)
   ns = numerical_setup(symbolic_setup(s,D_j),D_j)
 
-  b = GridapSolvers.allocate_col_vector(D_j)
-  x = GridapSolvers.allocate_col_vector(D_j)
-  y = GridapSolvers.allocate_row_vector(D_j)
+  b = GridapSolvers.allocate_in_domain(D_j)
+  x = GridapSolvers.allocate_in_domain(D_j)
+  y = GridapSolvers.allocate_in_range(D_j)
 
   fill!(b,1.0)
   y = b - D_j*x
@@ -24,8 +24,8 @@ end
 function test_solver(s,D_j)
   ns = numerical_setup(symbolic_setup(s,D_j),D_j)
 
-  b = GridapSolvers.allocate_col_vector(D_j)
-  x = GridapSolvers.allocate_col_vector(D_j)
+  b = GridapSolvers.allocate_in_domain(D_j)
+  x = GridapSolvers.allocate_in_domain(D_j)
 
   fill!(b,1.0)
   solve!(x,ns,b)
diff --git a/test/dev/Dj_solver_gmg.jl b/test/dev/Dj_solver_gmg.jl
index 3b36716..20155f5 100644
--- a/test/dev/Dj_solver_gmg.jl
+++ b/test/dev/Dj_solver_gmg.jl
@@ -105,8 +105,8 @@ end
 function test_solver(s,D_j)
   ns = numerical_setup(symbolic_setup(s,D_j),D_j)
 
-  b = GridapSolvers.allocate_col_vector(D_j)
-  x = GridapSolvers.allocate_col_vector(D_j)
+  b = GridapSolvers.allocate_in_domain(D_j)
+  x = GridapSolvers.allocate_in_domain(D_j)
 
   fill!(b,1.0)
   solve!(x,ns,b)
@@ -117,9 +117,9 @@ end
 
 function test_smoother(s,D_j)
   ns = numerical_setup(symbolic_setup(s,D_j),D_j)
-  b = GridapSolvers.allocate_col_vector(D_j)
-  x = GridapSolvers.allocate_col_vector(D_j)
-  r = GridapSolvers.allocate_row_vector(D_j)
+  b = GridapSolvers.allocate_in_domain(D_j)
+  x = GridapSolvers.allocate_in_domain(D_j)
+  r = GridapSolvers.allocate_in_range(D_j)
   fill!(b,1.0)
   fill!(x,1.0)
   mul!(r,D_j,x)
@@ -226,12 +226,12 @@ gmg = GMGLinearSolver(mh,
 solver = FGMRESSolver(100,gmg;rtol=1e-6,verbose=true)
 
 ns = numerical_setup(symbolic_setup(solver,A),A)
-x = GridapSolvers.allocate_col_vector(A)
+x = GridapSolvers.allocate_in_domain(A)
 solve!(x,ns,b)
 
 
 Pl = LinearSolvers.IdentitySolver()
 solver2 = GMRESSolver(1000;Pl=Pl,rtol=1e-6,verbose=true)
 ns2 = numerical_setup(symbolic_setup(solver2,A),A)
-x2 = GridapSolvers.allocate_col_vector(A)
+x2 = GridapSolvers.allocate_in_domain(A)
 solve!(x2,ns2,b)
diff --git a/test/dev/Dj_solver_petsc.jl b/test/dev/Dj_solver_petsc.jl
index 0c7bd4e..a168804 100644
--- a/test/dev/Dj_solver_petsc.jl
+++ b/test/dev/Dj_solver_petsc.jl
@@ -41,8 +41,8 @@ end
 function test_solver(s,D_j)
   ns = numerical_setup(symbolic_setup(s,D_j),D_j)
 
-  b = GridapSolvers.allocate_col_vector(D_j)
-  x = GridapSolvers.allocate_col_vector(D_j)
+  b = GridapSolvers.allocate_in_domain(D_j)
+  x = GridapSolvers.allocate_in_domain(D_j)
 
   fill!(b,1.0)
   solve!(x,ns,b)
diff --git a/test/dev/Lp_solver.jl b/test/dev/Lp_solver.jl
index af534f7..6338778 100644
--- a/test/dev/Lp_solver.jl
+++ b/test/dev/Lp_solver.jl
@@ -116,8 +116,8 @@ a_dg(p,v_p) = a_cg(p,v_p) + aΛ(p,v_p) + aΓ(p,v_p)
 Δp = assemble_matrix(a_dg,U_p,V_p)
 
 # Solver
-b = GridapSolvers.allocate_col_vector(Δp)
-x = GridapSolvers.allocate_col_vector(Δp)
+b = GridapSolvers.allocate_in_domain(Δp)
+x = GridapSolvers.allocate_in_domain(Δp)
 
 fill!(x,0.0)
 fill!(b,1.0)