From caffaeb2928e3ef1a52c567779aa266538402c4f Mon Sep 17 00:00:00 2001
From: Jesus Bonilla <jbonilla@cimne.upc.edu>
Date: Mon, 25 May 2020 17:47:51 +0200
Subject: [PATCH 1/3] Implemented periodic BC for CartesianModel (#263)

* Added field isperiodic to cartesian descriptor
* Modified constructors of CartesianDescriptor and CartesianDiscreteModel
* Added new method for _generate_cell_to_vertices_from_grid in
UnstructuredGrids.jl that implements the grid numbering for periodic BC.
* Added test PeriodicBC.jl in Geometry tests to check proper numbering of vefs
with periodic BC.
* Added tests PeriodicDarcy.jl and PeriodicCoupledPoisson.jl to check both
with periodic BC.
---
 src/Geometry/CartesianDiscreteModels.jl       |  10 +-
 src/Geometry/CartesianGrids.jl                | 128 ++++++++++++++--
 src/Geometry/UnstructuredGrids.jl             |  78 ++++++++++
 test/GeometryTests/PeriodicBCTests.jl         | 141 ++++++++++++++++++
 test/GeometryTests/runtests.jl                |   2 +
 .../PeriodicCoupledPoissonTests.jl            |  73 +++++++++
 test/GridapTests/PeriodicDarcyTests.jl        |  63 ++++++++
 test/GridapTests/runtests.jl                  |   4 +
 8 files changed, 484 insertions(+), 15 deletions(-)
 create mode 100644 test/GeometryTests/PeriodicBCTests.jl
 create mode 100644 test/GridapTests/PeriodicCoupledPoissonTests.jl
 create mode 100644 test/GridapTests/PeriodicDarcyTests.jl

diff --git a/src/Geometry/CartesianDiscreteModels.jl b/src/Geometry/CartesianDiscreteModels.jl
index 480861530..cb7db8457 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 = UnstructuredGridTopology(_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)
@@ -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..f2d537049 100644
--- a/src/Geometry/CartesianGrids.jl
+++ b/src/Geometry/CartesianGrids.jl
@@ -16,50 +16,154 @@ 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(
+      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,
+  isperiodic::NTuple{D,Bool}=tfill(false,Val{D}())) where D
+
+  T = eltype(sizes)
+  F = typeof(map)
+  CartesianDescriptor(origin,sizes,Tuple(partition),map,isperiodic)
+end
+
+"""
+    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,isperiodic)
 end
 
 """
     CartesianDescriptor(
-      pmin::Point{D},pmax::Point{D},partition,map::Function=identity) where D
+      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,
+  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,isperiodic)
+end
+
+"""
+    CartesianDescriptor(
+      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,
+  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,isperiodic)
+end
+
+"""
+    CartesianDescriptor(
+      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,isperiodic)
 end
 
 # Coordinates
@@ -181,12 +285,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
 
diff --git a/src/Geometry/UnstructuredGrids.jl b/src/Geometry/UnstructuredGrids.jl
index 25fc3cb3a..bd6daf131 100644
--- a/src/Geometry/UnstructuredGrids.jl
+++ b/src/Geometry/UnstructuredGrids.jl
@@ -167,6 +167,15 @@ function UnstructuredGridTopology(grid::UnstructuredGrid)
   _generate_grid_topology_from_grid(grid,cell_to_vertices,vertex_to_node)
 end
 
+function UnstructuredGridTopology(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 UnstructuredGridTopology(grid::UnstructuredGrid, cell_to_vertices::Table, vertex_to_node::AbstractVector)
   _generate_grid_topology_from_grid(grid,cell_to_vertices,vertex_to_node)
 end
@@ -213,6 +222,75 @@ function _generate_cell_to_vertices_from_grid(grid::UnstructuredGrid)
   (cell_to_vertices, vertex_to_node, node_to_vertex)
 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
+
+
 function _generate_cell_to_vertices(
   cell_to_nodes::Table,
   cell_to_cell_type::AbstractVector{<:Integer},
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

From 2fbbe2672990258e4d8d04c38847f7f8aed496f5 Mon Sep 17 00:00:00 2001
From: Jesus Bonilla <jbonilla@cimne.upc.edu>
Date: Tue, 26 May 2020 10:19:20 +0200
Subject: [PATCH 2/3] Updated CartesianDescriptor signatures (#263)

* Converted optional arguments of CartesianDescriptor constructors to
key-word arguments.
* Added deprecated signatures for backwards compatibility.
* Uptated tests to use new signatures.
* Updated News.md.
---
 NEWS.md                                       | 12 ++-
 src/Geometry/CartesianDiscreteModels.jl       |  2 +-
 src/Geometry/CartesianGrids.jl                | 83 +++----------------
 .../CartesianDiscreteModelsTests.jl           |  2 +-
 4 files changed, 22 insertions(+), 77 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index a4757f9c0..fa0ad6911 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,14 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [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).
 
 ## [0.10.2] - 2020-5-21
 
@@ -36,13 +44,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 cb7db8457..c8dbcce4f 100644
--- a/src/Geometry/CartesianDiscreteModels.jl
+++ b/src/Geometry/CartesianDiscreteModels.jl
@@ -47,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)
 
      grid = CartesianGrid(subdesc)
      _grid = UnstructuredGrid(grid)
diff --git a/src/Geometry/CartesianGrids.jl b/src/Geometry/CartesianGrids.jl
index f2d537049..1e0b6c17a 100644
--- a/src/Geometry/CartesianGrids.jl
+++ b/src/Geometry/CartesianGrids.jl
@@ -30,7 +30,7 @@ struct CartesianDescriptor{D,T,F<:Function} <: GridapType
   function CartesianDescriptor(
     origin::Point{D},
     sizes::NTuple{D},
-    partition,
+    partition;
     map::Function=identity,
     isperiodic::NTuple{D,Bool}=tfill(false,Val{D}())) where D
 
@@ -45,29 +45,6 @@ struct CartesianDescriptor{D,T,F<:Function} <: GridapType
     F = typeof(map)
     new{D,T,F}(origin,sizes,Tuple(partition),map,isperiodic)
   end
-
-end
-
-"""
-    CartesianDescriptor(
-      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,
-  isperiodic::NTuple{D,Bool}=tfill(false,Val{D}())) where D
-
-  T = eltype(sizes)
-  F = typeof(map)
-  CartesianDescriptor(origin,sizes,Tuple(partition),map,isperiodic)
 end
 
 """
@@ -77,28 +54,6 @@ end
       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,
-  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,isperiodic)
-end
-
-"""
-    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(
@@ -111,14 +66,14 @@ function CartesianDescriptor(
   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,isperiodic)
+  CartesianDescriptor(origin,sizes,partition;map=map,isperiodic=isperiodic)
 end
 
 """
     CartesianDescriptor(
       pmin::Point{D},
       pmax::Point{D},
-      partition,
+      partition;
       map::Function=identity,
       isperiodic::NTuple{D,Bool}=tfill(false,Val{D})) where D
 
@@ -127,7 +82,7 @@ end
 function CartesianDescriptor(
   pmin::Point{D},
   pmax::Point{D},
-  partition,
+  partition;
   map::Function=identity,
   isperiodic::NTuple{D,Bool}=tfill(false,Val{D}())) where D
 
@@ -137,34 +92,16 @@ function CartesianDescriptor(
     domain[2*(d-1)+1] = pmin[d]
     domain[2*(d-1)+2] = pmax[d]
   end
-  CartesianDescriptor(domain,partition,map,isperiodic)
+  CartesianDescriptor(domain,partition;map=map,isperiodic=isperiodic)
 end
 
-"""
-    CartesianDescriptor(
-      pmin::Point{D},
-      pmax::Point{D},
-      partition;
-      map::Function=identity,
-      isperiodic::NTuple{D,Bool}=tfill(false,Val{D}())) where D
+# Deprecated signatures for backwards compatibility
 
-`partition` is a 1D indexable collection of arbitrary type.
-"""
-function CartesianDescriptor(
-  pmin::Point{D},
-  pmax::Point{D},
-  partition;
-  map::Function=identity,
-  isperiodic::NTuple{D,Bool}=tfill(false,Val{D}())) where D
+@deprecate CartesianDescriptor(origin::Point{D}, sizes::NTuple{D}, partition, map::Function) where D CartesianDescriptor(origin, sizes, partition; map=map)
 
-  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,isperiodic)
-end
+@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
 
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
 

From 1fb958b6ba4e11c1df6994e6ad4f2ea3f248c621 Mon Sep 17 00:00:00 2001
From: Jesus Bonilla <jbonilla@cimne.upc.edu>
Date: Tue, 26 May 2020 15:19:42 +0200
Subject: [PATCH 3/3] Done changes suggested in PR #266

* Minor changes in CartesianDescriptor documentation
* Renamed UnstructuredGridTopology used for numbering a grid with periodic BC
to _cartesian_grid_topology_with_periodic_bcs.
* Moved above function (and related) to CartesianGrids.jl
---
 src/Geometry/CartesianDiscreteModels.jl |  4 +-
 src/Geometry/CartesianGrids.jl          | 83 ++++++++++++++++++++++++-
 src/Geometry/UnstructuredGrids.jl       | 80 +-----------------------
 3 files changed, 84 insertions(+), 83 deletions(-)

diff --git a/src/Geometry/CartesianDiscreteModels.jl b/src/Geometry/CartesianDiscreteModels.jl
index c8dbcce4f..cc0cda601 100644
--- a/src/Geometry/CartesianDiscreteModels.jl
+++ b/src/Geometry/CartesianDiscreteModels.jl
@@ -17,7 +17,7 @@ struct CartesianDiscreteModel{D,T,F} <: DiscreteModel{D,D}
     grid = CartesianGrid(desc)
     _grid = UnstructuredGrid(grid)
     if any(desc.isperiodic)
-      topo = UnstructuredGridTopology(_grid, desc.isperiodic, desc.partition)
+      topo = _cartesian_grid_topology_with_periodic_bcs(_grid, desc.isperiodic, desc.partition)
     else
       topo = UnstructuredGridTopology(_grid)
     end
@@ -47,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; map=desc.map)
+       CartesianDescriptor(Point(suborigin), subsizes, subpartition; map=desc.map, isperiodic=desc.isperiodic)
 
      grid = CartesianGrid(subdesc)
      _grid = UnstructuredGrid(grid)
diff --git a/src/Geometry/CartesianGrids.jl b/src/Geometry/CartesianGrids.jl
index 1e0b6c17a..a0f8663f0 100644
--- a/src/Geometry/CartesianGrids.jl
+++ b/src/Geometry/CartesianGrids.jl
@@ -21,7 +21,7 @@ struct CartesianDescriptor{D,T,F<:Function} <: GridapType
       CartesianDescriptor(
         origin::Point{D},
         sizes::NTuple{D},
-        partition,
+        partition;
         map::Function=identity,
         isperiodic::NTuple{D,Bool}=tfill(false,Val{D})) where D
 
@@ -50,7 +50,7 @@ end
 """
     CartesianDescriptor(
       domain,
-      partition,
+      partition;
       map::Function=identity,
       isperiodic::NTuple{D,Bool}=tfill(false,Val{D}))
 
@@ -272,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 bd6daf131..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)
@@ -167,15 +167,6 @@ function UnstructuredGridTopology(grid::UnstructuredGrid)
   _generate_grid_topology_from_grid(grid,cell_to_vertices,vertex_to_node)
 end
 
-function UnstructuredGridTopology(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 UnstructuredGridTopology(grid::UnstructuredGrid, cell_to_vertices::Table, vertex_to_node::AbstractVector)
   _generate_grid_topology_from_grid(grid,cell_to_vertices,vertex_to_node)
 end
@@ -222,74 +213,6 @@ function _generate_cell_to_vertices_from_grid(grid::UnstructuredGrid)
   (cell_to_vertices, vertex_to_node, node_to_vertex)
 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
-
 
 function _generate_cell_to_vertices(
   cell_to_nodes::Table,
@@ -419,4 +342,3 @@ function _refine_grid_connectivity(cell_to_points::Table, ltcell_to_lpoints)
     ltcell_to_lpoints)
   Table(data,ptrs)
 end
-