From d34e8a9df51644613f0bbd45b5f2b733023de6a8 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Thu, 15 Nov 2012 15:31:05 +0000
Subject: [PATCH] Mapping of Zaglmayr QuadBasis to fit Gmsh Ref Space

---
 FunctionSpace/QuadEdgeBasis.cpp | 27 +++++++++++++++++++++++++++
 FunctionSpace/QuadEdgeBasis.h   |  8 ++++++--
 FunctionSpace/QuadNodeBasis.cpp | 24 ++++++++++++++++++++++++
 FunctionSpace/QuadNodeBasis.h   |  8 ++++++--
 4 files changed, 63 insertions(+), 4 deletions(-)

diff --git a/FunctionSpace/QuadEdgeBasis.cpp b/FunctionSpace/QuadEdgeBasis.cpp
index 56e6570e3b..463e960fc1 100644
--- a/FunctionSpace/QuadEdgeBasis.cpp
+++ b/FunctionSpace/QuadEdgeBasis.cpp
@@ -201,6 +201,33 @@ QuadEdgeBasis::QuadEdgeBasis(int order){
   }
 
 
+  // Mapping to Gmsh Quad //
+  // x = (u + 1) / 2
+  // y = (v + 1) / 2
+  //
+  // (x, y) = Zaglmayr Ref Quad
+  // (u, v) = Gmsh     Ref Quad
+
+  Polynomial  mapX(Polynomial(0.5, 1, 0, 0) + 
+		   Polynomial(0.5, 0, 0, 0));
+
+  Polynomial  mapY(Polynomial(0.5, 0, 1, 0) + 
+		   Polynomial(0.5, 0, 0, 0));
+
+  for(int i = 0; i < nEdgeClosure; i++){
+    for(int j = 0; j < nEdge; j++){
+      (*(*edge)[i])[j]->at(0) = (*(*edge)[i])[j]->at(0).compose(mapX, mapY);
+      (*(*edge)[i])[j]->at(1) = (*(*edge)[i])[j]->at(1).compose(mapX, mapY);
+      (*(*edge)[i])[j]->at(2) = (*(*edge)[i])[j]->at(2).compose(mapX, mapY);
+    }
+  }
+
+  for(int i = 0; i < nCell; i++){
+    (*cell)[i]->at(0) = (*cell)[i]->at(0).compose(mapX, mapY);
+    (*cell)[i]->at(1) = (*cell)[i]->at(1).compose(mapX, mapY);
+    (*cell)[i]->at(2) = (*cell)[i]->at(2).compose(mapX, mapY);
+  }
+
   // Free Temporary Sapce //
   delete[] legendre;
   delete[] intLegendre;
diff --git a/FunctionSpace/QuadEdgeBasis.h b/FunctionSpace/QuadEdgeBasis.h
index a0400dd6de..40db8df987 100644
--- a/FunctionSpace/QuadEdgeBasis.h
+++ b/FunctionSpace/QuadEdgeBasis.h
@@ -10,10 +10,14 @@
    This class can instantiate an Edge-Based Basis 
    (high or low order) for Quads.@n
    
-   It uses 
+   It uses a variation of
    <a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a>  
    Basis for @em high @em order Polynomial%s generation.@n
- */
+
+   The following mapping has been applied to Zaglmayr's Basis for Quads:
+   @li @f$x = \frac{u + 1}{2}@f$
+   @li @f$y = \frac{v + 1}{2}@f$ 
+*/
 
 class QuadEdgeBasis: public BasisVector{
  public:
diff --git a/FunctionSpace/QuadNodeBasis.cpp b/FunctionSpace/QuadNodeBasis.cpp
index 68b0fc26a0..253d8f049e 100644
--- a/FunctionSpace/QuadNodeBasis.cpp
+++ b/FunctionSpace/QuadNodeBasis.cpp
@@ -127,6 +127,30 @@ QuadNodeBasis::QuadNodeBasis(int order){
   }
 
 
+  // Mapping to Gmsh Quad //
+  // x = (u + 1) / 2
+  // y = (v + 1) / 2
+  //
+  // (x, y) = Zaglmayr Ref Quad
+  // (u, v) = Gmsh     Ref Quad
+
+  Polynomial  mapX(Polynomial(0.5, 1, 0, 0) + 
+		   Polynomial(0.5, 0, 0, 0));
+
+  Polynomial  mapY(Polynomial(0.5, 0, 1, 0) + 
+		   Polynomial(0.5, 0, 0, 0));
+
+  for(int i = 0; i < nVertex; i++)
+    *(*node)[i] = (*node)[i]->compose(mapX, mapY);
+
+  for(int i = 0; i < nEdgeClosure; i++)
+    for(int j = 0; j < nEdge; j++)
+      *(*(*edge)[i])[j] = (*(*edge)[i])[j]->compose(mapX, mapY);
+
+  for(int i = 0; i < nCell; i++)
+    *(*cell)[i] = (*cell)[i]->compose(mapX, mapY);
+  
+
   // Free Temporary Sapce //
   delete[] legendre;
 }
diff --git a/FunctionSpace/QuadNodeBasis.h b/FunctionSpace/QuadNodeBasis.h
index 308aa78cf3..80f4ff0432 100644
--- a/FunctionSpace/QuadNodeBasis.h
+++ b/FunctionSpace/QuadNodeBasis.h
@@ -10,10 +10,14 @@
    This class can instantiate a Node-Based Basis 
    (high or low order) for Quads.@n
    
-   It uses 
+   It uses a variation of
    <a href="http://www.hpfem.jku.at/publications/szthesis.pdf">Zaglmayr's</a>  
    Basis for @em high @em order Polynomial%s generation.@n
- */
+ 
+   The following mapping has been applied to Zaglmayr's Basis for Quads:
+   @li @f$x = \frac{u + 1}{2}@f$
+   @li @f$y = \frac{v + 1}{2}@f$
+*/
 
 class QuadNodeBasis: public BasisScalar{
  public:
-- 
GitLab