diff --git a/fem/src/ElmerSolver.F90 b/fem/src/ElmerSolver.F90
index 1cf4cc98dd..8f710124ed 100644
--- a/fem/src/ElmerSolver.F90
+++ b/fem/src/ElmerSolver.F90
@@ -416,7 +416,10 @@ END SUBROUTINE TrilinosCleanup
          ! Optionally perform simple extrusion to increase the dimension of the mesh
          !----------------------------------------------------------------------------------
          CALL CreateExtrudedMesh() 
-
+         DO j=1,CurrentModel % NumberOfSolvers
+           CALL CreateExtrudedMesh(j) 
+         END DO
+         
          !----------------------------------------------------------------------------------
          ! If requested perform coordinate transformation directly after is has been obtained.
          ! Don't maintain the original mesh. 
@@ -843,30 +846,69 @@ END SUBROUTINE InitializeRandomSeed
      
      ! Optionally create extruded mesh on-the-fly.
      !--------------------------------------------------------------------
-     SUBROUTINE CreateExtrudedMesh()
-
+     SUBROUTINE CreateExtrudedMesh(SolverId)
+       INTEGER, OPTIONAL :: SolverId
+       
        LOGICAL :: SliceVersion
-
-       IF(.NOT. ListCheckPrefix(CurrentModel % Simulation,'Extruded Mesh') ) RETURN
+       TYPE(ValueList_t), POINTER :: VList
+       TYPE(Mesh_t), POINTER :: Mesh_in, pMesh, prevMesh
+       LOGICAL :: PrimaryMesh
+       
+       IF(PRESENT(SolverId)) THEN
+         Vlist => CurrentModel % Solvers(SolverId) % Values
+         Mesh_in => CurrentModel % Solvers(SolverId) % Mesh
+         PrimaryMesh = .FALSE.
+       ELSE
+         Vlist => CurrentModel % Simulation
+         Mesh_in => CurrentModel % Meshes
+         PrimaryMesh = .TRUE.
+       END IF
+                  
+       IF(.NOT. ListCheckPrefix(VList,'Extruded Mesh') ) RETURN
+       IF(.NOT. PrimaryMesh) THEN
+         CALL Info('CreateExtrudedMesh','Extruding mesh associated to solver '//I2S(SolverId))
+       END IF
        
-       SliceVersion = GetLogical(CurrentModel % Simulation,'Extruded Mesh Slices',Found )              
+       SliceVersion = GetLogical(Vlist,'Extruded Mesh Slices',Found )              
        IF( SliceVersion ) THEN
-         ExtrudedMesh => MeshExtrudeSlices(CurrentModel % Meshes, CurrentModel % Simulation )
+         ExtrudedMesh => MeshExtrudeSlices(Mesh_in, Vlist ) 
        ELSE
-         ExtrudedMesh => MeshExtrude(CurrentModel % Meshes, CurrentModel % Simulation)
+         ExtrudedMesh => MeshExtrude(Mesh_in, Vlist )
        END IF
          
        ! Make the solvers point to the extruded mesh, not the original mesh
        !-------------------------------------------------------------------
        DO i=1,CurrentModel % NumberOfSolvers
-         IF(ASSOCIATED(CurrentModel % Solvers(i) % Mesh,CurrentModel % Meshes)) &
-             CurrentModel % Solvers(i) % Mesh => ExtrudedMesh 
+         IF(ASSOCIATED(CurrentModel % Solvers(i) % Mesh,Mesh_in)) THEN
+           CALL Info('CreateExtrudedMesh','Pointing solver '//I2S(i)//' to the new extruded mesh!')
+           CurrentModel % Solvers(i) % Mesh => ExtrudedMesh
+         END IF
+       END DO
+
+       ! Put the extruded mesh in a correct place in the list of meshes.
+       i = 0
+       NULLIFY(prevMesh) 
+       pMesh => CurrentModel % Meshes
+       DO WHILE(ASSOCIATED(pMesh))
+         i = i+1
+         IF(ASSOCIATED(pMesh,Mesh_in)) EXIT
+         prevMesh => pMesh
+         pMesh => pMesh % Next 
        END DO
-       ExtrudedMesh % Next => CurrentModel % Meshes % Next
-       CurrentModel % Meshes => ExtrudedMesh
+       CALL Info('CreateExtrduedMesh','Extruded mesh order is '//I2S(i),Level=25)
 
+       ExtrudedMesh % Next => Mesh_in % Next
+       IF(ASSOCIATED(prevMesh)) THEN
+         prevMesh % Next => ExtrudedMesh
+       ELSE
+         CurrentModel % Meshes => ExtrudedMesh
+       END IF
+         
        ! If periodic BC given, compute boundary mesh projector:
+       ! but only for the primary mesh.
        ! ------------------------------------------------------
+       IF(.NOT. PrimaryMesh) RETURN
+       
        DO i = 1,CurrentModel % NumberOfBCs
          IF(ASSOCIATED(CurrentModel % Bcs(i) % PMatrix)) &
              CALL FreeMatrix( CurrentModel % BCs(i) % PMatrix )
diff --git a/fem/src/MeshUtils.F90 b/fem/src/MeshUtils.F90
index 9e9819a2e4..5746816a00 100644
--- a/fem/src/MeshUtils.F90
+++ b/fem/src/MeshUtils.F90
@@ -17900,7 +17900,7 @@ FUNCTION MeshExtrude(Mesh_in, Vlist) RESULT(Mesh_out)
     ! Create the division for the 1D mesh
     !--------------------------------------------
     CALL ExtrudedDivision(Vlist,nlev,Wtable)    
-    CALL Info(Caller,'Creating '//I2S(nlev)//' extruded element layers',Level=10)
+    CALL Info(Caller,'Extruding '//I2S(nlev)//' element layers on: '//TRIM(Mesh_in % Name),Level=10)
     in_levels = nlev-1
         
     ! Generate volume nodal points:
@@ -18602,7 +18602,8 @@ FUNCTION MeshExtrudeSlices(Mesh_in, Vlist) RESULT(Mesh_out)
     ! Create the division for the 1D mesh
     !--------------------------------------------
     CALL ExtrudedDivision(Vlist,nlev,Wtable)    
-    CALL Info(Caller,'Creating '//I2S(nlev)//' extruded element layers',Level=10)
+    CALL Info(Caller,'Extruding '//I2S(nlev)//' element layers on: '//TRIM(Mesh_in % Name),Level=10)
+
     
     ! In parallel let us pick only our own share of the
     ! division. This logic makes it possible to have nonuniform divisions easily.
diff --git a/fem/src/modules/Mesh2MeshSolver.F90 b/fem/src/modules/Mesh2MeshSolver.F90
index 893dbd3dd7..7bce1767ba 100644
--- a/fem/src/modules/Mesh2MeshSolver.F90
+++ b/fem/src/modules/Mesh2MeshSolver.F90
@@ -113,7 +113,8 @@ END SUBROUTINE InterpolateMeshToMesh
   CALL Info('Mesh2MeshSolver','Mapping result between meshes')
     
   ThisMesh => Getmesh()
-  CALL Info('Mesh2MeshSolver','This mesh name is: '//TRIM(ThisMesh % Name),Level=7)   
+  CALL Info('Mesh2MeshSolver','This mesh dimension is '&
+      //I2S(ThisMesh % MeshDim)//' and name: '//TRIM(ThisMesh % Name),Level=7)   
 
   Params => GetSolverParams()
 
@@ -136,6 +137,7 @@ END SUBROUTINE InterpolateMeshToMesh
         TargetMesh => Mesh
         EXIT
       END IF
+      Mesh => Mesh % Next
     END DO
     IF( ASSOCIATED( TargetMesh ) ) THEN
       CALL Info('Mesh2MeshSolver','Using target mesh as the first mesh different from this mesh',Level=8)