From 9ad2edf9af5a2e8a0216a6858d8526e1dc06e643 Mon Sep 17 00:00:00 2001
From: Nicolas Marsic <nicolas.marsic@gmail.com>
Date: Fri, 9 Nov 2012 13:55:11 +0000
Subject: [PATCH] Adaptive View for vectorial fields

---
 FunctionSpace/FunctionSpaceEdge.cpp | 37 ++++++++++++++++++++++
 FunctionSpace/FunctionSpaceEdge.h   |  5 +++
 FunctionSpace/FunctionSpaceNode.cpp | 20 ++++++++++++
 FunctionSpace/FunctionSpaceNode.h   |  5 +++
 FunctionSpace/FunctionSpaceScalar.h | 26 +++++++++++++++-
 FunctionSpace/FunctionSpaceVector.h | 26 +++++++++++++++-
 FunctionSpace/LagrangeBasis.cpp     | 48 ++++++++++++++++++++++-------
 FunctionSpace/LagrangeBasis.h       | 31 +++++++++++++++----
 8 files changed, 179 insertions(+), 19 deletions(-)

diff --git a/FunctionSpace/FunctionSpaceEdge.cpp b/FunctionSpace/FunctionSpaceEdge.cpp
index 1907ca1414..f7dc3b54ac 100644
--- a/FunctionSpace/FunctionSpaceEdge.cpp
+++ b/FunctionSpace/FunctionSpaceEdge.cpp
@@ -61,3 +61,40 @@ interpolate(const MElement& element,
   // Return Interpolated Value //
   return val;
 }
+
+fullVector<double> FunctionSpaceEdge::
+interpolateInRefSpace(const MElement& element, 
+		      const std::vector<double>& coef,
+		      const fullVector<double>& uvw) const{
+ 
+  // Const Cast For MElement //
+  MElement& eelement = 
+    const_cast<MElement&>(element);
+  
+  // Get Basis Functions //
+  const vector<const vector<Polynomial>*>& fun = getLocalFunctions(element);
+  const unsigned int nFun                      = fun.size();
+
+  // Get Jacobian //
+  fullMatrix<double>  invJac(3, 3);        
+  eelement.getJacobian(uvw(0), uvw(1), uvw(2), invJac);
+  invJac.invertInPlace();
+ 
+  // Interpolate (in Reference Place) //
+  fullVector<double> val(3); 
+  val(0) = 0; 
+  val(1) = 0; 
+  val(2) = 0;
+
+  for(unsigned int i = 0; i < nFun; i++){
+    fullVector<double> vi = 
+      Mapper::grad(Polynomial::at(*fun[i], uvw(0), uvw(1), uvw(2)),
+		   invJac);
+    
+    vi.scale(coef[i]);
+    val.axpy(vi, 1);
+  }
+
+  // Return Interpolated Value //
+  return val;
+}
diff --git a/FunctionSpace/FunctionSpaceEdge.h b/FunctionSpace/FunctionSpaceEdge.h
index bc03737d2c..0fbbaf9b11 100644
--- a/FunctionSpace/FunctionSpaceEdge.h
+++ b/FunctionSpace/FunctionSpaceEdge.h
@@ -21,6 +21,11 @@ class FunctionSpaceEdge : public FunctionSpaceVector{
     interpolate(const MElement& element, 
 		const std::vector<double>& coef,
 		const fullVector<double>& xyz) const;
+
+  virtual fullVector<double> 
+    interpolateInRefSpace(const MElement& element, 
+			  const std::vector<double>& coef,
+			  const fullVector<double>& uvw) const;
 };
 
 
diff --git a/FunctionSpace/FunctionSpaceNode.cpp b/FunctionSpace/FunctionSpaceNode.cpp
index 10928c769c..0adb1198b6 100644
--- a/FunctionSpace/FunctionSpaceNode.cpp
+++ b/FunctionSpace/FunctionSpaceNode.cpp
@@ -48,3 +48,23 @@ interpolate(const MElement& element,
   // Return Interpolated Value //
   return val;
 }
+
+double FunctionSpaceNode::
+interpolateInRefSpace(const MElement& element, 
+		      const std::vector<double>& coef,
+		      const fullVector<double>& uvw) const{
+
+  // Get Basis Functions //
+  const vector<const Polynomial*> fun = getLocalFunctions(element);
+  const unsigned int nFun             = fun.size();
+  
+  // Interpolate (in Reference Place) //
+  double val = 0; 
+
+  for(unsigned int i = 0; i < nFun; i++)
+    val += 
+      fun[i]->at(uvw(0), uvw(1), uvw(2)) * coef[i];
+
+  // Return Interpolated Value //
+  return val;
+}
diff --git a/FunctionSpace/FunctionSpaceNode.h b/FunctionSpace/FunctionSpaceNode.h
index 59655bcd48..d682c7da9e 100644
--- a/FunctionSpace/FunctionSpaceNode.h
+++ b/FunctionSpace/FunctionSpaceNode.h
@@ -21,6 +21,11 @@ class FunctionSpaceNode : public FunctionSpaceScalar{
     interpolate(const MElement& element, 
 		const std::vector<double>& coef,
 		const fullVector<double>& xyz) const;
+
+  virtual double 
+    interpolateInRefSpace(const MElement& element, 
+			  const std::vector<double>& coef,
+			  const fullVector<double>& uvw) const;
 };
 
 
diff --git a/FunctionSpace/FunctionSpaceScalar.h b/FunctionSpace/FunctionSpaceScalar.h
index 18a86c9773..c51c27458b 100644
--- a/FunctionSpace/FunctionSpaceScalar.h
+++ b/FunctionSpace/FunctionSpaceScalar.h
@@ -35,6 +35,11 @@ class FunctionSpaceScalar : public FunctionSpace{
     interpolate(const MElement& element, 
 		const std::vector<double>& coef,
 		const fullVector<double>& xyz) const = 0;
+
+  virtual double 
+    interpolateInRefSpace(const MElement& element, 
+			  const std::vector<double>& coef,
+			  const fullVector<double>& uvw) const = 0;
   
   const std::vector<const Polynomial*> 
     getLocalFunctions(const MElement& element) const;
@@ -57,7 +62,26 @@ class FunctionSpaceScalar : public FunctionSpace{
    @param coef The coefficients of the interpolation
    @param xyz The coordinate 
    (of point @em inside the given @c element)
-   of the interpolation
+   of the interpolation in the @em Physical Space
+
+   @return Returns the (scalar) interpolated value
+
+   @warning
+   If the given coordinate are not in the given
+   @c element @em Bad @em Things may happend
+   
+   @todo
+   If the given coordinate are not in the given
+   @c element @em Bad @em Things may happend
+   ---> check
+   **
+
+   @fn FunctionSpaceScalar::interpolateInRefSpace
+   @param element The MElement to interpolate on
+   @param coef The coefficients of the interpolation
+   @param uvw The coordinate 
+   (of point @em inside the given @c element)
+   of the interpolation in the @em Reference Space
 
    @return Returns the (scalar) interpolated value
 
diff --git a/FunctionSpace/FunctionSpaceVector.h b/FunctionSpace/FunctionSpaceVector.h
index 3b4ff01119..7227f306b2 100644
--- a/FunctionSpace/FunctionSpaceVector.h
+++ b/FunctionSpace/FunctionSpaceVector.h
@@ -40,6 +40,11 @@ class FunctionSpaceVector : public FunctionSpace{
 		const std::vector<double>& coef,
 		const fullVector<double>& xyz) const = 0;
 
+  virtual fullVector<double> 
+    interpolateInRefSpace(const MElement& element, 
+			  const std::vector<double>& coef,
+			  const fullVector<double>& uvw) const = 0;
+
   const std::vector<const std::vector<Polynomial>*>
     getLocalFunctions(const MElement& element) const;
 
@@ -64,7 +69,26 @@ class FunctionSpaceVector : public FunctionSpace{
    @param coef The coefficients of the interpolation
    @param xyz The coordinate 
    (of point @em inside the given @c element)
-   of the interpolation
+   of the interpolation in the @em Physical Space
+
+   @return Returns the (vectorial) interpolated value
+
+   @warning
+   If the given coordinate are not in the given
+   @c element @em Bad @em Things may happend
+   
+   @todo
+   If the given coordinate are not in the given
+   @c element @em Bad @em Things may happend
+   ---> check
+   **
+
+   @fn FunctionSpaceVector::interpolateInRefSpace
+   @param element The MElement to interpolate on
+   @param coef The coefficients of the interpolation
+   @param xyz The coordinate 
+   (of point @em inside the given @c element)
+   of the interpolation in the @em Reference Space
 
    @return Returns the (vectorial) interpolated value
 
diff --git a/FunctionSpace/LagrangeBasis.cpp b/FunctionSpace/LagrangeBasis.cpp
index 63b51e1efe..9402df5f4c 100644
--- a/FunctionSpace/LagrangeBasis.cpp
+++ b/FunctionSpace/LagrangeBasis.cpp
@@ -9,24 +9,50 @@ LagrangeBasis::LagrangeBasis(void){
 LagrangeBasis::~LagrangeBasis(void){
 }
 
-fullVector<double> LagrangeBasis::
-project(const fullVector<double>& coef,
-	const std::vector<const Polynomial*>& basis){
+vector<double> LagrangeBasis::
+project(const MElement& element,
+	const vector<double>& coef,
+	const FunctionSpaceScalar& fSpace){
 
   // Init New Coefs //
   const unsigned int size = point->size1();
-  fullVector<double> newCoef(size);
+  vector<double> newCoef(size);
 
   // Interpolation at Lagrange Points //
-  const unsigned int nCoef = coef.size();
- 
   for(unsigned int i = 0; i < size; i++){
-    newCoef(i) = 0;
+    fullVector<double> uvw(3);
+    uvw(0) = (*point)(i, 0);
+    uvw(1) = (*point)(i, 1);
+    uvw(2) = (*point)(i, 2);
+
+    newCoef[i] = fSpace.interpolateInRefSpace(element, 
+					      coef, 
+					      uvw);
+  }
+    
+  // Return ;
+  return newCoef;
+}
+
+vector<fullVector<double> > LagrangeBasis::
+project(const MElement& element,
+	const vector<double>& coef,
+	const FunctionSpaceVector& fSpace){
+  
+  // Init New Coefs //
+  const unsigned int size = point->size1();
+  vector<fullVector<double> > newCoef(size);
+
+  // Interpolation at Lagrange Points //
+  for(unsigned int i = 0; i < size; i++){
+    fullVector<double> uvw(3);
+    uvw(0) = (*point)(i, 0);
+    uvw(1) = (*point)(i, 1);
+    uvw(2) = (*point)(i, 2);
     
-    for(unsigned int j = 0; j < nCoef; j++)
-      newCoef(i) += coef(j) * basis[j]->at((*point)(i, 0),
-					   (*point)(i, 1),
-					   (*point)(i, 2));
+    newCoef[i] = fSpace.interpolateInRefSpace(element, 
+					      coef, 
+					      uvw);
   }
     
   // Return ;
diff --git a/FunctionSpace/LagrangeBasis.h b/FunctionSpace/LagrangeBasis.h
index 0d29cb6c78..7daa7f49fe 100644
--- a/FunctionSpace/LagrangeBasis.h
+++ b/FunctionSpace/LagrangeBasis.h
@@ -6,6 +6,10 @@
 #include "polynomialBasis.h"
 #include "BasisScalar.h"
 
+#include "MElement.h"
+#include "FunctionSpaceScalar.h"
+#include "FunctionSpaceVector.h"
+
 /**
    @interface LagrangeBasis
    @brief Interoface for Lagrange Basis
@@ -39,13 +43,28 @@ class LagrangeBasis: public BasisScalar{
   //! @return Returns the Monomial Matrix
   const fullMatrix<double>& getMonomial(void) const;
 
-  //! @param coef A vector of Real numbers
-  //! @param basis A vector of Polynomials 
-  //! in a @em Reference Space
+  //! @param element A MElement
+  //! @param coef A vector of coefficient associated 
+  //! to the given Element
+  //! @param fSpace The (scalar) Function Space 
+  //! of the given Coefficients
   //! @return Projects the given Coefficients in this LagrangeBasis@n
-  fullVector<double> project(const fullVector<double>& coef,
-			     const std::vector<const Polynomial*>& basis);
-
+  std::vector<double> project(const MElement& element,
+			      const std::vector<double>& coef,
+			      const FunctionSpaceScalar& fSpace);
+
+  //! @param element A MElement
+  //! @param coef A vector of coefficient associated 
+  //! to the given Element
+  //! @param fSpace The (vectorial) Function Space 
+  //! of the given Coefficients
+  //! @return Projects the given Coefficients in this LagrangeBasis@n
+  //! @note Each Coefficients will be projected into a vector
+  //! with the same dimesion as the vectorial Polynomials
+  std::vector<fullVector<double> > 
+    project(const MElement& element,
+	    const std::vector<double>& coef,
+	    const FunctionSpaceVector& fSpace);
  protected:
   //! Returns a new LagrangeBasis
   //!
-- 
GitLab