diff --git a/NEWS.md b/NEWS.md
index dde54a996..333e5016f 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -6,6 +6,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+ - Extended support of `CartesianDiscreteModel` to models with periodic boundary conditions.
+ PR [#266](https://github.com/gridap/Gridap.jl/pull/266).
+
+### Deprecated
+ - Optional argument `map` for CartesianDescriptor converted to a key-word argument. Since PR [#266](https://github.com/gridap/Gridap.jl/pull/266).
+ 
 ### Fixed
  - Fixed some methods of the `sparsecsr` generic function. Since PR [#262](https://github.com/gridap/Gridap.jl/pull/262).
  - Fixed BUG in `findnz` function for `SparseMatrixCSR`. Since PR [#264](https://github.com/gridap/Gridap.jl/pull/264).
@@ -41,13 +48,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   - The part associated with the imposition of constraints in the `FESpace` interface has changed slightly. Since PR [#247](https://github.com/gridap/Gridap.jl/pull/247).
   - Simplified the signature of `zero_free_values(::FESpace)`. Since PR [#249](https://github.com/gridap/Gridap.jl/pull/249).
   - Simplified the signature of `zero_initial_guess(op::NonlinearOperator)`. Since PR [#249](https://github.com/gridap/Gridap.jl/pull/249).
-  - Major refactoring in the `Assembler` interface. 
+  - Major refactoring in the `Assembler` interface.
     **Important change:** Now, assembly-related functions take the data returned by functions like
     `collect_cell_matrix` as it is. Example: the old user code `assemble_matrix(assembler,collect_cell_matrix(du,dv,terms)...)`
     now is written simply as `assemble_matrix(assembler,collect_cell_matrix(du,dv,terms))`, i.e., the unpack of the last argument is not
     used anymore.  In addition, with the new assembler interface, it is possible to customize the assembly process
     via a so-called `AssemblerStrategy` object. Since PR [#249](https://github.com/gridap/Gridap.jl/pull/249).
-  - Change the types of the sizes and partition fields of CartesianDescriptor to tuples instead of points. 
+  - Change the types of the sizes and partition fields of CartesianDescriptor to tuples instead of points.
     Since PR [#246](https://github.com/gridap/Gridap.jl/pull/246).
 
 ## [0.9.2] - 2020-4-26
diff --git a/src/Geometry/CartesianDiscreteModels.jl b/src/Geometry/CartesianDiscreteModels.jl
index 480861530..cc0cda601 100644
--- a/src/Geometry/CartesianDiscreteModels.jl
+++ b/src/Geometry/CartesianDiscreteModels.jl
@@ -16,7 +16,11 @@ struct CartesianDiscreteModel{D,T,F} <: DiscreteModel{D,D}
   function CartesianDiscreteModel(desc::CartesianDescriptor{D,T,F}) where {D,T,F}
     grid = CartesianGrid(desc)
     _grid = UnstructuredGrid(grid)
-    topo = UnstructuredGridTopology(_grid)
+    if any(desc.isperiodic)
+      topo = _cartesian_grid_topology_with_periodic_bcs(_grid, desc.isperiodic, desc.partition)
+    else
+      topo = UnstructuredGridTopology(_grid)
+    end
     nfaces = [num_faces(topo,d) for d in 0:num_cell_dims(topo)]
     labels = FaceLabeling(nfaces)
     _fill_cartesian_face_labeling!(labels,topo)
@@ -43,7 +47,7 @@ struct CartesianDiscreteModel{D,T,F} <: DiscreteModel{D,D}
      subpartition = Tuple(cmax) .- Tuple(cmin) .+ 1
      subsizes = desc.sizes
      subdesc =
-       CartesianDescriptor(Point(suborigin), subsizes, subpartition, desc.map)
+       CartesianDescriptor(Point(suborigin), subsizes, subpartition; map=desc.map, isperiodic=desc.isperiodic)
 
      grid = CartesianGrid(subdesc)
      _grid = UnstructuredGrid(grid)
@@ -60,8 +64,8 @@ end
 
 Same args needed to construct a `CartesianDescriptor`
 """
-function CartesianDiscreteModel(args...)
-  desc = CartesianDescriptor(args...)
+function CartesianDiscreteModel(args...; kwargs...)
+  desc = CartesianDescriptor(args...; kwargs...)
   CartesianDiscreteModel(desc)
 end
 
diff --git a/src/Geometry/CartesianGrids.jl b/src/Geometry/CartesianGrids.jl
index 356694fbb..a0f8663f0 100644
--- a/src/Geometry/CartesianGrids.jl
+++ b/src/Geometry/CartesianGrids.jl
@@ -16,52 +16,93 @@ struct CartesianDescriptor{D,T,F<:Function} <: GridapType
   sizes::NTuple{D,T}
   partition::NTuple{D,Int}
   map::F
+  isperiodic::NTuple{D,Bool}
   @doc """
       CartesianDescriptor(
-        origin::Point{D}, sizes::NTuple{D}, partition, map::Function=identity) where D
+        origin::Point{D},
+        sizes::NTuple{D},
+        partition;
+        map::Function=identity,
+        isperiodic::NTuple{D,Bool}=tfill(false,Val{D})) where D
 
   `partition` is a 1D indexable collection of arbitrary type.
   """
   function CartesianDescriptor(
-    origin::Point{D}, sizes::NTuple{D}, partition, map::Function=identity) where D
+    origin::Point{D},
+    sizes::NTuple{D},
+    partition;
+    map::Function=identity,
+    isperiodic::NTuple{D,Bool}=tfill(false,Val{D}())) where D
+
+    for i in 1:D
+      if isperiodic[i]
+        @assert (partition[i] > 2) "A minimum of 3 elements is required in any
+          periodic direction"
+      end
+    end
 
     T = eltype(sizes)
     F = typeof(map)
-    new{D,T,F}(origin,sizes,Tuple(partition),map)
+    new{D,T,F}(origin,sizes,Tuple(partition),map,isperiodic)
   end
-
 end
 
 """
-    CartesianDescriptor(domain,partition,map::Function=identity)
+    CartesianDescriptor(
+      domain,
+      partition;
+      map::Function=identity,
+      isperiodic::NTuple{D,Bool}=tfill(false,Val{D}))
 
 `domain` and `partition` are 1D indexable collections of arbitrary type.
 """
-function CartesianDescriptor(domain,partition,map::Function=identity)
+function CartesianDescriptor(
+  domain,
+  partition;
+  map::Function=identity,
+  isperiodic::NTuple=tfill(false,Val{length(partition)}()))
+
   D = length(partition)
   limits = [(domain[2*d-1],domain[2*d]) for d in 1:D]
   sizes = Tuple([(limits[d][2]-limits[d][1])/partition[d] for d in 1:D])
   origin = Point([ limits[d][1] for d in 1:D]...)
-  CartesianDescriptor(origin,sizes,partition,map)
+  CartesianDescriptor(origin,sizes,partition;map=map,isperiodic=isperiodic)
 end
 
 """
     CartesianDescriptor(
-      pmin::Point{D},pmax::Point{D},partition,map::Function=identity) where D
+      pmin::Point{D},
+      pmax::Point{D},
+      partition;
+      map::Function=identity,
+      isperiodic::NTuple{D,Bool}=tfill(false,Val{D})) where D
 
 `partition` is a 1D indexable collection of arbitrary type.
 """
 function CartesianDescriptor(
-  pmin::Point{D},pmax::Point{D},partition,map::Function=identity) where D
+  pmin::Point{D},
+  pmax::Point{D},
+  partition;
+  map::Function=identity,
+  isperiodic::NTuple{D,Bool}=tfill(false,Val{D}())) where D
+
   T = eltype(pmin)
   domain = zeros(T,2*D)
   for d in 1:D
     domain[2*(d-1)+1] = pmin[d]
     domain[2*(d-1)+2] = pmax[d]
   end
-  CartesianDescriptor(domain,partition,map)
+  CartesianDescriptor(domain,partition;map=map,isperiodic=isperiodic)
 end
 
+# Deprecated signatures for backwards compatibility
+
+@deprecate CartesianDescriptor(origin::Point{D}, sizes::NTuple{D}, partition, map::Function) where D CartesianDescriptor(origin, sizes, partition; map=map)
+
+@deprecate CartesianDescriptor(domain,partition,map::Function) CartesianDescriptor(domain,partition;map=map)
+
+@deprecate CartesianDescriptor(pmin::Point{D}, pmax::Point{D}, partition, map::Function) where D CartesianDescriptor(pmin, pmax, partition; map=map)
+
 # Coordinates
 
 struct CartesianCoordinates{D,T,F} <: AbstractArray{Point{D,T},D}
@@ -181,12 +222,12 @@ get_reffes(g::CartesianGrid{2}) = [QUAD4,]
 get_reffes(g::CartesianGrid{3}) = [HEX8,]
 
 """
-    CartesianGrid(args...)
+    CartesianGrid(args...;kwargs...)
 
 Same args needed to construct a `CartesianDescriptor`
 """
-function CartesianGrid(args...)
-  desc = CartesianDescriptor(args...)
+function CartesianGrid(args...;kwargs...)
+  desc = CartesianDescriptor(args...;kwargs...)
   CartesianGrid(desc)
 end
 
@@ -231,3 +272,82 @@ end
 function get_cell_map(grid::CartesianGrid{D,T,typeof(identity)} where {D,T})
   CartesianMap(grid.node_coords.data)
 end
+
+# Cartesian grid topology with periodic BC
+
+function _cartesian_grid_topology_with_periodic_bcs(grid::UnstructuredGrid,
+  isperiodic::NTuple,
+  partition)
+
+  cell_to_vertices, vertex_to_node =
+    _generate_cell_to_vertices_from_grid(grid, isperiodic, partition)
+  _generate_grid_topology_from_grid(grid,cell_to_vertices,vertex_to_node)
+end
+
+function _generate_cell_to_vertices_from_grid(grid::UnstructuredGrid,
+  isperiodic::NTuple, partition)
+
+  if is_first_order(grid)
+    nodes = get_cell_nodes(grid)
+    cell_to_vertices = copy(nodes)
+
+    nnodes = num_nodes(grid)
+    num_nodes_x_dir = [partition[i]+1 for i in 1:length(partition)]
+    point_to_isperiodic, slave_point_to_point, slave_point_to_master_point =
+      _generate_slave_to_master_point(num_nodes_x_dir,isperiodic, nnodes)
+
+    vertex_to_point = findall( .! point_to_isperiodic)
+    point_to_vertex = fill(-1,length(point_to_isperiodic))
+    point_to_vertex[vertex_to_point] = 1:length(vertex_to_point)
+    point_to_vertex[slave_point_to_point] = point_to_vertex[slave_point_to_master_point]
+
+    cell_to_vertices = Table(LocalToGlobalArray(nodes,point_to_vertex))
+
+    vertex_to_node = vertex_to_point
+    node_to_vertex = point_to_vertex
+  else
+    @notimplemented
+  end
+  (cell_to_vertices,vertex_to_node, node_to_vertex)
+end
+
+function _generate_slave_to_master_point(num_nodes_x_dir::Vector{Int},
+  isperiodic::NTuple, num_nodes::Int)
+
+  point_to_isperiodic = fill(false,num_nodes)
+
+  slave_ijk_bounds = Array{Any,1}(undef,length(isperiodic))
+  master_ijk_bounds = Array{Any,1}(undef,length(isperiodic))
+
+  linear_indices = LinearIndices(Tuple(num_nodes_x_dir))
+  periodic_dirs = findall(x->x==true, isperiodic)
+  for periodic_dir in periodic_dirs
+    for i in 1:length(isperiodic)
+      if i == periodic_dir
+        slave_ijk_bounds[i] = num_nodes_x_dir[i]
+      else
+        slave_ijk_bounds[i] = 1:num_nodes_x_dir[i]
+      end
+    end
+    slave_ijk_set = CartesianIndices(Tuple(slave_ijk_bounds))
+    point_to_isperiodic[linear_indices[slave_ijk_set]] .= true
+  end
+
+  slave_point_to_point = findall( point_to_isperiodic)
+  slave_point_to_master_point = Array{Int,1}(undef,length(slave_point_to_point))
+
+  cartesian_indices = CartesianIndices(Tuple(num_nodes_x_dir))
+  for (i,point) in enumerate(slave_point_to_point)
+    ijk = collect(cartesian_indices[point].I)
+    for i in periodic_dirs
+      if ijk[i] == num_nodes_x_dir[i]
+        ijk[i] = 1
+      end
+    end
+
+    master_point_ijk = CartesianIndex(Tuple(ijk))
+    slave_point_to_master_point[i] = linear_indices[master_point_ijk]
+  end
+
+  point_to_isperiodic, slave_point_to_point, slave_point_to_master_point
+end
diff --git a/src/Geometry/UnstructuredGrids.jl b/src/Geometry/UnstructuredGrids.jl
index 25fc3cb3a..459b8bbda 100644
--- a/src/Geometry/UnstructuredGrids.jl
+++ b/src/Geometry/UnstructuredGrids.jl
@@ -33,7 +33,7 @@ struct UnstructuredGrid{Dc,Dp,Tp,O} <: Grid{Dc,Dp}
   end
 end
 
-""" 
+"""
     UnstructuredGrid(trian::Grid)
 """
 function UnstructuredGrid(trian::Grid)
@@ -213,6 +213,7 @@ function _generate_cell_to_vertices_from_grid(grid::UnstructuredGrid)
   (cell_to_vertices, vertex_to_node, node_to_vertex)
 end
 
+
 function _generate_cell_to_vertices(
   cell_to_nodes::Table,
   cell_to_cell_type::AbstractVector{<:Integer},
@@ -341,4 +342,3 @@ function _refine_grid_connectivity(cell_to_points::Table, ltcell_to_lpoints)
     ltcell_to_lpoints)
   Table(data,ptrs)
 end
-
diff --git a/test/GeometryTests/CartesianDiscreteModelsTests.jl b/test/GeometryTests/CartesianDiscreteModelsTests.jl
index 8728a48bf..55639b9b7 100644
--- a/test/GeometryTests/CartesianDiscreteModelsTests.jl
+++ b/test/GeometryTests/CartesianDiscreteModelsTests.jl
@@ -36,7 +36,7 @@ end
 
 domain = (1,2,0,pi,0,0.5)
 partition = (10,30,4)
-model = CartesianDiscreteModel(domain,partition,polar)
+model = CartesianDiscreteModel(domain,partition;map=polar)
 test_discrete_model(model)
 @test is_oriented(get_grid(model)) == true
 
diff --git a/test/GeometryTests/PeriodicBCTests.jl b/test/GeometryTests/PeriodicBCTests.jl
new file mode 100644
index 000000000..f876f15d9
--- /dev/null
+++ b/test/GeometryTests/PeriodicBCTests.jl
@@ -0,0 +1,141 @@
+module PeriodicBCTests
+
+using Gridap.Arrays
+using Gridap.Geometry
+
+using Test
+
+
+## 2D mesh with periodic BC in x
+model = CartesianDiscreteModel((0,1,0,1),(3,3);isperiodic=(true,false))
+
+# Test cells to vertices
+@test model.grid_topology.n_m_to_nface_to_mfaces[3,1] ==
+  [[1, 2, 4, 5],
+   [2, 3, 5, 6],
+   [3, 1, 6, 4],
+   [4, 5, 7, 8],
+   [5, 6, 8, 9],
+   [6, 4, 9, 7],
+   [7, 8, 10, 11],
+   [8, 9, 11, 12],
+   [9, 7, 12, 10]]
+
+# Test cells to faces
+@test model.grid_topology.n_m_to_nface_to_mfaces[3,2] ==
+  [[1, 2, 3, 4],
+   [5, 6, 4, 7],
+   [8, 9, 7, 3],
+   [2, 10, 11, 12],
+   [6, 13, 12, 14],
+   [9, 15, 14, 11],
+   [10, 16, 17, 18],
+   [13, 19, 18, 20],
+   [15, 21, 20, 17]]
+
+# Test faces to vertices
+@test model.grid_topology.n_m_to_nface_to_mfaces[2,1] ==
+  [[1, 2], [4, 5], [1, 4], [2, 5], [2, 3], [5, 6], [3, 6], [3, 1], [6, 4],
+   [7, 8], [4, 7], [5, 8], [8, 9], [6, 9], [9, 7],
+   [10, 11], [7, 10], [8, 11], [11, 12], [9, 12], [12, 10]]
+
+# Test vertices to faces
+@test model.grid_topology.n_m_to_nface_to_mfaces[1,2] ==
+  [[1, 3, 8],
+   [1, 4, 5],
+   [5, 7, 8],
+   [2, 3, 9, 11],
+   [2, 4, 6, 12],
+   [6, 7, 9, 14],
+   [10, 11, 15, 17],
+   [10, 12, 13, 18],
+   [13, 14, 15, 20],
+   [16, 17, 21],
+   [16, 18, 19],
+   [19, 20, 21]]
+
+# Test vertices to cells
+@test model.grid_topology.n_m_to_nface_to_mfaces[1,3] ==
+  [[1, 3],
+   [1, 2],
+   [2, 3],
+   [1, 3, 4, 6],
+   [1, 2, 4, 5],
+   [2, 3, 5, 6],
+   [4, 6, 7, 9],
+   [4, 5, 7, 8],
+   [5, 6, 8, 9],
+   [7, 9],
+   [7, 8],
+   [8, 9]]
+
+# Test faces to cells
+@test model.grid_topology.n_m_to_nface_to_mfaces[2,3] ==
+  [[1], [1, 4], [1, 3], [1, 2], [2], [2, 5], [2, 3], [3], [3, 6],
+   [4, 7], [4, 6], [4, 5], [5, 8], [5, 6], [6, 9],
+   [7], [7, 9], [7, 8], [8], [8, 9], [9]]
+
+## 2D mesh with periodic BC in x & y
+model = CartesianDiscreteModel((0,1,0,1),(3,3);isperiodic=(true,true))
+
+ # Test cells to vertices
+@test model.grid_topology.n_m_to_nface_to_mfaces[3,1] ==
+  [[1, 2, 4, 5],
+   [2, 3, 5, 6],
+   [3, 1, 6, 4],
+   [4, 5, 7, 8],
+   [5, 6, 8, 9],
+   [6, 4, 9, 7],
+   [7, 8, 1, 2],
+   [8, 9, 2, 3],
+   [9, 7, 3, 1]]
+
+# Test cells to faces
+@test model.grid_topology.n_m_to_nface_to_mfaces[3,2] ==
+  [[1, 2, 3, 4],
+   [5, 6, 4, 7],
+   [8, 9, 7, 3],
+   [2, 10, 11, 12],
+   [6, 13, 12, 14],
+   [9, 15, 14, 11],
+   [10, 1, 16, 17],
+   [13, 5, 17, 18],
+   [15, 8, 18, 16]]
+
+# Test faces to vertices
+@test model.grid_topology.n_m_to_nface_to_mfaces[2,1] ==
+  [[1, 2], [4, 5], [1, 4], [2, 5], [2, 3], [5, 6], [3, 6], [3, 1], [6, 4],
+   [7, 8], [4, 7], [5, 8], [8, 9], [6, 9], [9, 7],
+   [7, 1], [8, 2], [9, 3]]
+
+# Test vertices to faces
+@test model.grid_topology.n_m_to_nface_to_mfaces[1,2] ==
+  [[1, 3, 8, 16],
+   [1, 4, 5, 17],
+   [5, 7, 8, 18],
+   [2, 3, 9, 11],
+   [2, 4, 6, 12],
+   [6, 7, 9, 14],
+   [10, 11, 15, 16],
+   [10, 12, 13, 17],
+   [13, 14, 15, 18]]
+
+# Test vertices to cells
+@test model.grid_topology.n_m_to_nface_to_mfaces[1,3] ==
+  [[1, 3, 7, 9],
+   [1, 2, 7, 8],
+   [2, 3, 8, 9],
+   [1, 3, 4, 6],
+   [1, 2, 4, 5],
+   [2, 3, 5, 6],
+   [4, 6, 7, 9],
+   [4, 5, 7, 8],
+   [5, 6, 8, 9]]
+
+# Test faces to cells
+@test model.grid_topology.n_m_to_nface_to_mfaces[2,3] ==
+  [[1, 7], [1, 4], [1, 3], [1, 2], [2, 8], [2, 5], [2, 3], [3, 9], [3, 6],
+   [4, 7], [4, 6], [4, 5], [5, 8], [5, 6], [6, 9],
+   [7, 9], [7, 8], [8, 9]]
+
+end #module
diff --git a/test/GeometryTests/runtests.jl b/test/GeometryTests/runtests.jl
index 448856d2e..1f38de550 100644
--- a/test/GeometryTests/runtests.jl
+++ b/test/GeometryTests/runtests.jl
@@ -36,6 +36,8 @@ using Test
 
 @testset "CartesianDiscreteModels" begin include("CartesianDiscreteModelsTests.jl") end
 
+@testset "PeriodicBC" begin include("PeriodicBCTests.jl") end
+
 @testset "BoundaryTriangulations" begin include("BoundaryTriangulationsTests.jl") end
 
 @testset "GenericBoundaryTriangulations" begin include("GenericBoundaryTriangulationsTests.jl") end
diff --git a/test/GridapTests/PeriodicCoupledPoissonTests.jl b/test/GridapTests/PeriodicCoupledPoissonTests.jl
new file mode 100644
index 000000000..bf6a97af3
--- /dev/null
+++ b/test/GridapTests/PeriodicCoupledPoissonTests.jl
@@ -0,0 +1,73 @@
+module PeriodicCoupledPoissonTests
+
+using Gridap
+using Test
+
+
+u(x) =  x[1]^2 + 2*x[2]^2
+v(x) = -x[2]^2
+f(x) = -Δ(u)(x)
+g(x) = -Δ(v)(x) - u(x)
+
+model = CartesianDiscreteModel((0,1,0,1,0,0.01),(4,4,3);isperiodic=(false,false,true))
+order = 2
+labels = get_face_labeling(model)
+
+add_tag_from_tags!(labels,"dirichlet",append!(collect(1:20),[23,24,25,26]))
+
+trian = get_triangulation(model)
+degree = 2*order
+quad = CellQuadrature(trian,degree)
+
+Vu = FESpace(
+       reffe=:Lagrangian, order=order, valuetype=Float64,
+       conformity=:H1, model=model, dirichlet_tags="dirichlet")
+Vv = FESpace(
+      reffe=:Lagrangian, order=order, valuetype=Float64,
+      conformity=:H1, model=model, dirichlet_tags="dirichlet")
+
+U = TrialFESpace(Vu,u)
+V = TrialFESpace(Vv,v)
+
+X = MultiFieldFESpace([U, V])
+Y = MultiFieldFESpace([Vu, Vv])
+
+function aa(X,Y)
+  u, v = X
+  v_u, v_v = Y
+  inner(∇(v_u),∇(u)) + inner(∇(v_v),∇(v)) - v_v*u
+end
+
+function l(Y)
+  v_u, v_v = Y
+  v_u*f + v_v*g
+end
+
+t_Ω = AffineFETerm(aa,l,trian,quad)
+op = AffineFEOperator(X,Y,t_Ω)
+uh, vh = solve(op)
+
+
+eu = u - uh
+ev = v - vh
+
+l2(u) = u*u
+h1(u) = ∇(u)*∇(u) + l2(u)
+
+eul2 = sqrt(sum( integrate(l2(eu),trian,quad) ))
+euh1 = sqrt(sum( integrate(h1(eu),trian,quad) ))
+ul2 = sqrt(sum( integrate(l2(uh),trian,quad) ))
+uh1 = sqrt(sum( integrate(h1(uh),trian,quad) ))
+
+evl2 = sqrt(sum( integrate(l2(ev),trian,quad) ))
+evh1 = sqrt(sum( integrate(h1(ev),trian,quad) ))
+vl2 = sqrt(sum( integrate(l2(vh),trian,quad) ))
+vh1 = sqrt(sum( integrate(h1(vh),trian,quad) ))
+
+@test eul2/ul2 < 1.e-8
+@test euh1/uh1 < 1.e-7
+@test evl2/vl2 < 1.e-8
+@test evh1/vh1 < 1.e-7
+
+
+end #module
diff --git a/test/GridapTests/PeriodicDarcyTests.jl b/test/GridapTests/PeriodicDarcyTests.jl
new file mode 100644
index 000000000..65c4d4787
--- /dev/null
+++ b/test/GridapTests/PeriodicDarcyTests.jl
@@ -0,0 +1,63 @@
+module PeriodicDarcyTests
+
+using Gridap
+using Test
+
+u(x) = VectorValue(x[1]*(x[1]-1)*(2x[2]-1.0),-x[2]*(x[2]-1.0)*(2x[1]-1.0))
+p(x) = x[2]-0.5
+f(x) = VectorValue(x[1]*(x[1]-1)*(2x[2]-1.0),-x[2]*(x[2]-1.0)*(2x[1]-1.0)+1.0)
+g(x) = 0.0
+
+domain = (0,1,0,1,0)
+order = 1
+partition = (4,4)
+model = CartesianDiscreteModel(domain,partition; isperiodic=(true,false))
+
+V = FESpace(
+reffe=:RaviartThomas, order=order, valuetype=VectorValue{2,Float64},
+conformity=:Hdiv, model=model,dirichlet_tags=collect(1:6))
+
+Q = FESpace(
+reffe=:QLagrangian, order=order, valuetype=Float64,
+conformity=:L2, model=model, constraint=:zeromean)
+
+U = TrialFESpace(V,u)
+P = TrialFESpace(Q)
+
+Y = MultiFieldFESpace([V, Q])
+X = MultiFieldFESpace([U, P])
+
+trian = Triangulation(model)
+degree = 2*(order)
+quad = CellQuadrature(trian,degree)
+x = get_physical_coordinate(trian)
+
+function a(x,y)
+  u, p = x
+  v, q = y
+  v*u - p*(∇*v) + q*(∇*u)
+end
+
+function l(y)
+  v, q = y
+  v*f + q*g
+end
+
+t_Ω = AffineFETerm(a,l,trian,quad)
+op = AffineFEOperator(X,Y,t_Ω)
+xh = solve(op)
+uh, ph = xh
+
+eu = u - uh
+ep = p - ph
+
+l2(v) = v*v
+eu_l2 = sum(integrate(l2(eu),trian,quad))
+ep_l2 = sum(integrate(l2(ep),trian,quad))
+
+tol = 1.0e-9
+@test eu_l2 < tol
+@test ep_l2 < tol
+
+
+end # module
diff --git a/test/GridapTests/runtests.jl b/test/GridapTests/runtests.jl
index c9e11e633..384fc6dab 100644
--- a/test/GridapTests/runtests.jl
+++ b/test/GridapTests/runtests.jl
@@ -16,6 +16,10 @@ using Test
 
 @testset "Darcy" begin include("DarcyTests.jl") end
 
+@testset "PeriodicDarcy" begin include("PeriodicDarcyTests.jl") end
+
+@testset "PeriodicCoupledPoisson" begin include("PeriodicCoupledPoissonTests.jl") end
+
 @testset "SurfaceCoupling" begin include("SurfaceCouplingTests.jl") end
 
 @testset "IsotropicDamage" begin include("IsotropicDamageTests.jl") end