From 8b6d39e386b7221b8c88f68e3733e977afff660a Mon Sep 17 00:00:00 2001
From: Francesc Verdugo <fverdugo@cimne.upc.edu>
Date: Tue, 17 Mar 2020 10:04:49 +0100
Subject: [PATCH] Changed the definition of order in Raviart-Thomas and Nedelec
 reffes

---
 src/Polynomials/QCurlGradMonomialBases.jl     | 10 +++-
 src/Polynomials/QGradMonomialBases.jl         | 28 +++++++----
 src/ReferenceFEs/NedelecRefFEs.jl             | 15 +++---
 src/ReferenceFEs/RaviartThomasRefFEs.jl       | 11 +++--
 .../CurlConformingFESpacesTests.jl            |  2 +-
 .../DivConformingFESpacesTests.jl             |  2 +-
 test/GridapTests/DarcyTests.jl                |  6 +--
 .../QCurlGradMonomialBasesTests.jl            | 19 +++++++-
 .../QGradMonomialBasesTests.jl                | 20 +++++++-
 test/ReferenceFEsTests/NedelecRefFEsTests.jl  | 46 +++++++++++++++++++
 .../RaviartThomasRefFEsTests.jl               | 17 ++++++-
 test/ReferenceFEsTests/runtests.jl            |  2 +
 12 files changed, 146 insertions(+), 32 deletions(-)
 create mode 100644 test/ReferenceFEsTests/NedelecRefFEsTests.jl

diff --git a/src/Polynomials/QCurlGradMonomialBases.jl b/src/Polynomials/QCurlGradMonomialBases.jl
index e042a4ee3..222ac90fb 100644
--- a/src/Polynomials/QCurlGradMonomialBases.jl
+++ b/src/Polynomials/QCurlGradMonomialBases.jl
@@ -21,11 +21,14 @@ end
 
 Returns a `QCurlGradMonomialBasis` object. `D` is the dimension
 of the coordinate space and `T` is the type of the components in the vector-value.
+The `order` argument has the following meaning: the divergence of the  functions in this basis
+is in the Q space of degree `order`.
 """
 function QCurlGradMonomialBasis{D}(::Type{T},order::Int) where {D,T}
   @assert T<:Real "T needs to be <:Real since represents the type of the components of the vector value"
-  _t = tfill(order,Val{D-1}())
-  t = (order+1,_t...)
+  _order = order+1
+  _t = tfill(_order,Val{D-1}())
+  t = (_order+1,_t...)
   terms = CartesianIndices(t)
   perms = _prepare_perms(D)
   QCurlGradMonomialBasis(T,order,terms,perms)
@@ -53,3 +56,6 @@ get_value_type(::QCurlGradMonomialBasis{D,T}) where {D,T} = T
     num_terms(f::QCurlGradMonomialBasis{D,T}) where {D,T}
 """
 num_terms(f::QCurlGradMonomialBasis{D,T}) where {D,T} = length(f.qgrad.terms)*D
+
+get_order(f::QCurlGradMonomialBasis{D,T}) where {D,T} = get_order(f.qgrad)
+
diff --git a/src/Polynomials/QGradMonomialBases.jl b/src/Polynomials/QGradMonomialBases.jl
index 972fb693d..c2139a4e7 100644
--- a/src/Polynomials/QGradMonomialBases.jl
+++ b/src/Polynomials/QGradMonomialBases.jl
@@ -22,21 +22,31 @@ end
 
 Returns a `QGradMonomialBasis` object. `D` is the dimension
 of the coordinate space and `T` is the type of the components in the vector-value.
+The `order` argument has the following meaning: the curl of the  functions in this basis
+is in the Q space of degree `order`.
 """
 function QGradMonomialBasis{D}(::Type{T},order::Int) where {D,T}
   @assert T<:Real "T needs to be <:Real since represents the type of the components of the vector value"
-  _t = tfill(order+1,Val{D-1}())
-  t = (order,_t...)
+  _order = order + 1
+  _t = tfill(_order+1,Val{D-1}())
+  t = (_order,_t...)
   terms = CartesianIndices(t)
   perms = _prepare_perms(D)
   QGradMonomialBasis(T,order,terms,perms)
 end
 
+"""
+    num_terms(f::QGradMonomialBasis{D,T}) where {D,T}
+"""
+num_terms(f::QGradMonomialBasis{D,T}) where {D,T} = length(f.terms)*D
+
+get_order(f::QGradMonomialBasis) = f.order
+
 function field_cache(f::QGradMonomialBasis{D,T},x) where {D,T}
   @assert D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
   ndof = _ndofs_qgrad(f)
-  n = 1 + f.order
+  n = 1 + f.order+1
   V = VectorValue{D,T}
   r = CachedArray(zeros(V,(np,ndof)))
   v = CachedArray(zeros(V,(ndof,)))
@@ -48,13 +58,13 @@ function evaluate_field!(cache,f::QGradMonomialBasis{D,T},x) where {D,T}
   r, v, c = cache
   np = length(x)
   ndof = _ndofs_qgrad(f)
-  n = 1 + f.order
+  n = 1 + f.order+1
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
   setsize!(c,(D,n))
   for i in 1:np
     @inbounds xi = x[i]
-      _evaluate_nd_qgrad!(v,xi,f.order,f.terms,f.perms,c)
+      _evaluate_nd_qgrad!(v,xi,f.order+1,f.terms,f.perms,c)
     for j in 1:ndof
       @inbounds r[i,j] = v[j]
     end
@@ -66,7 +76,7 @@ function gradient_cache(f::QGradMonomialBasis{D,T},x) where {D,T}
   @assert D == length(eltype(x)) "Incorrect number of point components"
   np = length(x)
   ndof = _ndofs_qgrad(f)
-  n = 1 + f.order
+  n = 1 + f.order+1
   xi = testitem(x)
   V = VectorValue{D,T}
   G = gradient_type(V,xi)
@@ -81,7 +91,7 @@ function evaluate_gradient!(cache,f::QGradMonomialBasis{D,T},x) where {D,T}
   r, v, c, g = cache
   np = length(x)
   ndof = _ndofs_qgrad(f)
-  n = 1 + f.order
+  n = 1 + f.order+1
   setsize!(r,(np,ndof))
   setsize!(v,(ndof,))
   setsize!(c,(D,n))
@@ -89,7 +99,7 @@ function evaluate_gradient!(cache,f::QGradMonomialBasis{D,T},x) where {D,T}
   V = VectorValue{D,T}
   for i in 1:np
     @inbounds xi = x[i]
-    _gradient_nd_qgrad!(v,xi,f.order,f.terms,f.perms,c,g,V)
+    _gradient_nd_qgrad!(v,xi,f.order+1,f.terms,f.perms,c,g,V)
     for j in 1:ndof
       @inbounds r[i,j] = v[j]
     end
@@ -99,8 +109,6 @@ end
 
 # Helpers
 
-#_ndofs_qgrad(f::QGradMonomialBasis{D}) where D = D*f.order*(f.order+1)^(D-1)
-
 _ndofs_qgrad(f::QGradMonomialBasis{D}) where D = D*(length(f.terms))
 
 function _prepare_perms(D)
diff --git a/src/ReferenceFEs/NedelecRefFEs.jl b/src/ReferenceFEs/NedelecRefFEs.jl
index c0161f08c..25b8c5f38 100644
--- a/src/ReferenceFEs/NedelecRefFEs.jl
+++ b/src/ReferenceFEs/NedelecRefFEs.jl
@@ -1,5 +1,8 @@
 """
     NedelecRefFE(::Type{et},p::Polytope,order::Integer) where et
+
+The `order` argument has the following meaning: the curl of the  functions in this basis
+is in the Q space of degree `order`.
 """
 # @santiagobadia : Project, go to complex numbers
 function NedelecRefFE(::Type{et},p::Polytope,order::Integer) where et
@@ -48,7 +51,7 @@ function _Nedelec_nodes_and_moments(::Type{et}, p::Polytope, order::Integer) whe
   nf_nodes[erange] = ecips
   nf_moments[erange] = emoments
 
-  if ( num_dims(p) == 3 && order > 1)
+  if ( num_dims(p) == 3 && order > 0)
 
     fcips, fmoments = _Nedelec_face_values(p,et,order)
     frange = get_dimrange(p,D-1)
@@ -57,7 +60,7 @@ function _Nedelec_nodes_and_moments(::Type{et}, p::Polytope, order::Integer) whe
 
   end
 
-  if (order > 1)
+  if (order > 0)
 
     ccips, cmoments = _Nedelec_cell_values(p,et,order)
     crange = get_dimrange(p,D)
@@ -79,7 +82,7 @@ function _Nedelec_edge_values(p,et,order)
   egeomap = _ref_face_to_faces_geomap(p,ep)
 
   # Compute integration points at all polynomial edges
-  degree = order*2
+  degree = (order+1)*2
   equad = Quadrature(ep,degree)
   cips = get_coordinates(equad)
   wips = get_weights(equad)
@@ -88,7 +91,7 @@ function _Nedelec_edge_values(p,et,order)
   c_eips, ecips, ewips = _nfaces_evaluation_points_weights(p, egeomap, cips, wips)
 
   # Edge moments, i.e., M(Ei)_{ab} = q_RE^a(xgp_REi^b) w_Fi^b t_Ei ⋅ ()
-  eshfs = MonomialBasis(et,ep,order-1)
+  eshfs = MonomialBasis(et,ep,order)
   emoments = _Nedelec_edge_moments(p, eshfs, c_eips, ecips, ewips)
 
   return ecips, emoments
@@ -116,7 +119,7 @@ function _Nedelec_face_values(p,et,order)
   fgeomap = _ref_face_to_faces_geomap(p,fp)
 
   # Compute integration points at all polynomial edges
-  degree = order*2
+  degree = (order+1)*2
   fquad = Quadrature(fp,degree)
   fips = get_coordinates(fquad)
   wips = get_weights(fquad)
@@ -156,7 +159,7 @@ end
 function _Nedelec_cell_values(p,et,order)
 
   # Compute integration points at interior
-  degree = 2*order
+  degree = 2*(order+1)
   iquad = Quadrature(p,degree)
   ccips = get_coordinates(iquad)
   cwips = get_weights(iquad)
diff --git a/src/ReferenceFEs/RaviartThomasRefFEs.jl b/src/ReferenceFEs/RaviartThomasRefFEs.jl
index db4c4f9ba..87fb7f8c8 100644
--- a/src/ReferenceFEs/RaviartThomasRefFEs.jl
+++ b/src/ReferenceFEs/RaviartThomasRefFEs.jl
@@ -1,6 +1,9 @@
 
 """
     RaviartThomasRefFE(::Type{et},p::Polytope,order::Integer) where et
+
+The `order` argument has the following meaning: the divergence of the  functions in this basis
+is in the Q space of degree `order`.
 """
 function RaviartThomasRefFE(::Type{et},p::Polytope,order::Integer) where et
 
@@ -48,7 +51,7 @@ function _RT_nodes_and_moments(::Type{et}, p::Polytope, order::Integer) where et
   nf_nodes[frange] = fcips
   nf_moments[frange] = fmoments
 
-  if (order > 1)
+  if (order > 0)
     ccips, cmoments = _RT_cell_values(p,et,order)
     crange = get_dimrange(p,D)
     nf_nodes[crange] = ccips
@@ -115,7 +118,7 @@ function _RT_face_values(p,et,order)
   # the one of the points in the polytope after applying the geopmap
   # (fcips), and the weights for these nodes (fwips, a constant cell array)
   # Nodes (fcips)
-  degree = order*2
+  degree = (order+1)*2
   fquad = Quadrature(fp,degree)
   fips = get_coordinates(fquad)
   wips = get_weights(fquad)
@@ -124,7 +127,7 @@ function _RT_face_values(p,et,order)
 
   # Moments (fmoments)
   # The RT prebasis is expressed in terms of shape function
-  fshfs = MonomialBasis(et,fp,order-1)
+  fshfs = MonomialBasis(et,fp,order)
 
   # Face moments, i.e., M(Fi)_{ab} = q_RF^a(xgp_RFi^b) w_Fi^b n_Fi ⋅ ()
   fmoments = _RT_face_moments(p, fshfs, c_fips, fcips, fwips)
@@ -142,7 +145,7 @@ end
 # It provides for every cell the nodes and the moments arrays
 function _RT_cell_values(p,et,order)
   # Compute integration points at interior
-  degree = 2*order
+  degree = 2*(order+1)
   iquad = Quadrature(p,degree)
   ccips = get_coordinates(iquad)
   cwips = get_weights(iquad)
diff --git a/test/FESpacesTests/CurlConformingFESpacesTests.jl b/test/FESpacesTests/CurlConformingFESpacesTests.jl
index 8a87ead71..9b451018e 100644
--- a/test/FESpacesTests/CurlConformingFESpacesTests.jl
+++ b/test/FESpacesTests/CurlConformingFESpacesTests.jl
@@ -10,7 +10,7 @@ partition = (3,3,3)
 model = CartesianDiscreteModel(domain,partition)
 trian = get_triangulation(model)
 
-order = 3
+order = 2
 
 u(x) = VectorValue(x[1]*x[1],x[1]*x[1]*x[1],0.0)
 # u(x) = x
diff --git a/test/FESpacesTests/DivConformingFESpacesTests.jl b/test/FESpacesTests/DivConformingFESpacesTests.jl
index 34a9fe2bb..4c8e5e437 100644
--- a/test/FESpacesTests/DivConformingFESpacesTests.jl
+++ b/test/FESpacesTests/DivConformingFESpacesTests.jl
@@ -11,7 +11,7 @@ partition = (3,3)
 model = CartesianDiscreteModel(domain,partition)
 trian = get_triangulation(model)
 
-order = 2
+order = 1
 
 u(x) = x
 
diff --git a/test/GridapTests/DarcyTests.jl b/test/GridapTests/DarcyTests.jl
index 70362667e..97b620370 100644
--- a/test/GridapTests/DarcyTests.jl
+++ b/test/GridapTests/DarcyTests.jl
@@ -18,7 +18,7 @@ f(x) = u(x) + ∇p(x)
 
 domain = (0,1,0,1)
 partition = (4,4)
-order = 2
+order = 1
 model = CartesianDiscreteModel(domain,partition)
 
 V = FESpace(
@@ -26,7 +26,7 @@ V = FESpace(
   conformity=:Hdiv, model=model, dirichlet_tags=[5,6])
 
 Q = FESpace(
-  reffe=:QLagrangian, order=order-1, valuetype=Float64,
+  reffe=:QLagrangian, order=order, valuetype=Float64,
   conformity=:L2, model=model)
 
 U = TrialFESpace(V,u)
@@ -41,7 +41,7 @@ quad = CellQuadrature(trian,degree)
 
 neumanntags = [7,8]
 btrian = BoundaryTriangulation(model,neumanntags)
-degree = 2*order
+degree = 2*(order+1)
 bquad = CellQuadrature(btrian,degree)
 nb = get_normal_vector(btrian)
 
diff --git a/test/PolynomialsTests/QCurlGradMonomialBasesTests.jl b/test/PolynomialsTests/QCurlGradMonomialBasesTests.jl
index 5b95558f6..5dbb6c812 100644
--- a/test/PolynomialsTests/QCurlGradMonomialBasesTests.jl
+++ b/test/PolynomialsTests/QCurlGradMonomialBasesTests.jl
@@ -1,14 +1,29 @@
 module QCurlGradMonomialBasesTests
 
+using Test
 using Gridap.TensorValues
 using Gridap.Fields
 using Gridap.Polynomials
 
+xi = Point(2,3)
+np = 5
+x = fill(xi,np)
+
+order = 0
+D = 2
+T = Float64
+V = VectorValue{D,T}
+G = gradient_type(V,xi)
+b = QCurlGradMonomialBasis{D}(T,order)
+
+@test num_terms(b) == 4
+@test get_order(b) == 0
+
 xi = Point(2,3,5)
 np = 5
 x = fill(xi,np)
 
-order = 1
+order = 0
 D = 3
 T = Float64
 V = VectorValue{D,T}
@@ -35,7 +50,7 @@ xi = Point(2,3)
 np = 5
 x = fill(xi,np)
 
-order = 2
+order = 1
 D = 2
 T = Float64
 V = VectorValue{D,T}
diff --git a/test/PolynomialsTests/QGradMonomialBasesTests.jl b/test/PolynomialsTests/QGradMonomialBasesTests.jl
index d99bf1b44..47db170a8 100644
--- a/test/PolynomialsTests/QGradMonomialBasesTests.jl
+++ b/test/PolynomialsTests/QGradMonomialBasesTests.jl
@@ -1,14 +1,30 @@
 module QGradMonomialBasesTests
 
+using Test
 using Gridap.TensorValues
 using Gridap.Fields
 using Gridap.Polynomials
 
+xi = Point(2,3)
+np = 5
+x = fill(xi,np)
+
+order = 0
+D = 2
+T = Float64
+V = VectorValue{D,T}
+G = gradient_type(V,xi)
+b = QGradMonomialBasis{D}(T,order)
+
+@test num_terms(b) == 4
+@test b.order == 0
+@test get_order(b) == 0
+
 xi = Point(2,3,5)
 np = 5
 x = fill(xi,np)
 
-order = 1
+order = 0
 D = 3
 T = Float64
 V = VectorValue{D,T}
@@ -43,7 +59,7 @@ xi = Point(2,3)
 np = 5
 x = fill(xi,np)
 
-order = 2
+order = 1
 D = 2
 T = Float64
 V = VectorValue{D,T}
diff --git a/test/ReferenceFEsTests/NedelecRefFEsTests.jl b/test/ReferenceFEsTests/NedelecRefFEsTests.jl
new file mode 100644
index 000000000..20e356022
--- /dev/null
+++ b/test/ReferenceFEsTests/NedelecRefFEsTests.jl
@@ -0,0 +1,46 @@
+module NedelecRefFEsTest
+
+using Test
+using Gridap.Polynomials
+using Gridap.Fields
+using Gridap.TensorValues
+using Gridap.Fields: MockField
+using Gridap.ReferenceFEs
+
+p = QUAD
+D = num_dims(QUAD)
+et = Float64
+order = 0
+
+reffe = NedelecRefFE(et,p,order)
+test_reference_fe(reffe)
+@test num_terms(get_prebasis(reffe)) == 4
+@test get_order(get_prebasis(reffe)) == 0
+@test num_dofs(reffe) == 4
+
+p = QUAD
+D = num_dims(QUAD)
+et = Float64
+order = 1
+
+reffe = NedelecRefFE(et,p,order)
+test_reference_fe(reffe)
+@test num_terms(get_prebasis(reffe)) == 12
+@test num_dofs(reffe) == 12
+@test get_order(get_prebasis(reffe)) == 1
+
+prebasis = get_prebasis(reffe)
+dof_basis = get_dof_basis(reffe)
+
+v = VectorValue(3.0,0.0)
+field = MockField{D}(v)
+
+cache = dof_cache(dof_basis,field)
+r = evaluate_dof!(cache, dof_basis, field)
+test_dof(dof_basis,field,r)
+
+cache = dof_cache(dof_basis,prebasis)
+r = evaluate_dof!(cache, dof_basis, prebasis)
+test_dof(dof_basis,prebasis,r)
+
+end # module
diff --git a/test/ReferenceFEsTests/RaviartThomasRefFEsTests.jl b/test/ReferenceFEsTests/RaviartThomasRefFEsTests.jl
index b3bac0304..64e2bfe3c 100644
--- a/test/ReferenceFEsTests/RaviartThomasRefFEsTests.jl
+++ b/test/ReferenceFEsTests/RaviartThomasRefFEsTests.jl
@@ -1,5 +1,6 @@
 module RaviartThomasRefFEsTest
 
+using Test
 using Gridap.Polynomials
 using Gridap.Fields
 using Gridap.TensorValues
@@ -9,10 +10,24 @@ using Gridap.ReferenceFEs
 p = QUAD
 D = num_dims(QUAD)
 et = Float64
-order = 3
+order = 0
 
 reffe = RaviartThomasRefFE(et,p,order)
 test_reference_fe(reffe)
+@test num_terms(get_prebasis(reffe)) == 4
+@test get_order(get_prebasis(reffe)) == 0
+@test num_dofs(reffe) == 4
+
+p = QUAD
+D = num_dims(QUAD)
+et = Float64
+order = 1
+
+reffe = RaviartThomasRefFE(et,p,order)
+test_reference_fe(reffe)
+@test num_terms(get_prebasis(reffe)) == 12
+@test num_dofs(reffe) == 12
+@test get_order(get_prebasis(reffe)) == 1
 
 prebasis = get_prebasis(reffe)
 dof_basis = get_dof_basis(reffe)
diff --git a/test/ReferenceFEsTests/runtests.jl b/test/ReferenceFEsTests/runtests.jl
index 385af6b2f..34752e82e 100644
--- a/test/ReferenceFEsTests/runtests.jl
+++ b/test/ReferenceFEsTests/runtests.jl
@@ -24,4 +24,6 @@ using Test
 
 @testset "RaviartThomasRefFEs" begin include("RaviartThomasRefFEsTests.jl") end
 
+@testset "NedelecRefFEs" begin include("NedelecRefFEsTests.jl") end
+
 end # module