Skip to content
Snippets Groups Projects
Commit 2c89bb93 authored by Nicolas Marsic's avatar Nicolas Marsic
Browse files

Cleaning -- Works for Tris -- Reference Space *Error* for *Quads*

parent eddff855
No related branches found
No related tags found
No related merge requests found
...@@ -25,7 +25,7 @@ class BasisScalar: public Basis{ ...@@ -25,7 +25,7 @@ class BasisScalar: public Basis{
//! @return Returns the set of @em Polynomial%s //! @return Returns the set of @em Polynomial%s
//! defining this (scalar) Basis //! defining this (scalar) Basis
const std::vector<Polynomial>& getBasis(void) const; const std::vector<Polynomial>& getFunctions(void) const;
protected: protected:
//! Instantiate a new BasisScalar //! Instantiate a new BasisScalar
...@@ -37,7 +37,10 @@ class BasisScalar: public Basis{ ...@@ -37,7 +37,10 @@ class BasisScalar: public Basis{
// Inline Functions // // Inline Functions //
////////////////////// //////////////////////
inline const std::vector<Polynomial>& BasisScalar::getBasis(void) const{ inline
const std::vector<Polynomial>& BasisScalar::
getFunctions(void) const{
return *basis; return *basis;
} }
......
...@@ -25,7 +25,7 @@ class BasisVector: public Basis{ ...@@ -25,7 +25,7 @@ class BasisVector: public Basis{
//! @return Returns the set of @em Polynomial%s //! @return Returns the set of @em Polynomial%s
//! defining this (vectorial) Basis //! defining this (vectorial) Basis
const std::vector<std::vector<Polynomial> >& getBasis(void) const; const std::vector<std::vector<Polynomial> >& getFunctions(void) const;
protected: protected:
//! Instantiate a new BasisVector //! Instantiate a new BasisVector
...@@ -37,7 +37,10 @@ class BasisVector: public Basis{ ...@@ -37,7 +37,10 @@ class BasisVector: public Basis{
// Inline Functions // // Inline Functions //
////////////////////// //////////////////////
inline const std::vector<std::vector<Polynomial> >& BasisVector::getBasis(void) const{ inline
const std::vector<std::vector<Polynomial> >& BasisVector::
getFunctions(void) const{
return *basis; return *basis;
} }
......
...@@ -40,7 +40,6 @@ class FunctionSpace{ ...@@ -40,7 +40,6 @@ class FunctionSpace{
virtual ~FunctionSpace(void); virtual ~FunctionSpace(void);
const GroupOfElement& getSupport(void) const; const GroupOfElement& getSupport(void) const;
const Basis& getBasis(const MElement& element) const;
int getType(void) const; int getType(void) const;
unsigned int getNFunctionPerVertex(const MElement& element) const; unsigned int getNFunctionPerVertex(const MElement& element) const;
...@@ -71,10 +70,6 @@ inline const GroupOfElement& FunctionSpace::getSupport(void) const{ ...@@ -71,10 +70,6 @@ inline const GroupOfElement& FunctionSpace::getSupport(void) const{
return *goe; return *goe;
} }
inline const Basis& FunctionSpace::getBasis(const MElement& element) const{
return *basis;
}
inline int FunctionSpace::getType(void) const{ inline int FunctionSpace::getType(void) const{
return type; return type;
} }
......
...@@ -15,38 +15,38 @@ FunctionSpaceEdge::FunctionSpaceEdge(const GroupOfElement& goe, ...@@ -15,38 +15,38 @@ FunctionSpaceEdge::FunctionSpaceEdge(const GroupOfElement& goe,
FunctionSpaceEdge::~FunctionSpaceEdge(void){ FunctionSpaceEdge::~FunctionSpaceEdge(void){
} }
#include <cstdio>
fullVector<double> FunctionSpaceEdge:: fullVector<double> FunctionSpaceEdge::
interpolate(const MElement& element, interpolate(const MElement& element,
const std::vector<double>& coef, const std::vector<double>& coef,
const fullVector<double>& xyz) const{ const fullVector<double>& xyz) const{
// Const Cast For MElement //
MElement& eelement = MElement& eelement =
const_cast<MElement&>(element); const_cast<MElement&>(element);
const vector<vector<Polynomial> >& fun = // Get Basis Functions //
static_cast<const BasisVector&>(*basis).getBasis(); const vector<vector<Polynomial> >& fun = getBasis(element).getFunctions();
const unsigned int nFun = fun.size(); const unsigned int nFun = fun.size();
fullVector<double> val(3); // Get Jacobian //
fullMatrix<double> invJac(3, 3); fullMatrix<double> invJac(3, 3);
eelement.getJacobian(xyz(0), xyz(1), xyz(2), invJac); eelement.getJacobian(xyz(0), xyz(1), xyz(2), invJac);
invJac.invertInPlace(); invJac.invertInPlace();
// Get First Vertex of Element //
fullVector<double> origin(3); fullVector<double> origin(3);
origin(0) = eelement.getVertex(0)->x(); origin(0) = eelement.getVertex(0)->x();
origin(1) = eelement.getVertex(0)->y(); origin(1) = eelement.getVertex(0)->y();
origin(2) = eelement.getVertex(0)->z(); origin(2) = eelement.getVertex(0)->z();
// Map XYZ to Reference Plane (UVW) //
fullVector<double> uvw = Mapper::invMap(xyz, origin, invJac); fullVector<double> uvw = Mapper::invMap(xyz, origin, invJac);
printf("%g, %g, %g\n", uvw(0), uvw(1), uvw(2)); // Interpolate (in Reference Place) //
fullVector<double> val(3);
//val(0) = 1; val(1) = 1; val(2) = 0; val(0) = 0;
val(1) = 0;
val(2) = 0;
for(unsigned int i = 0; i < nFun; i++){ for(unsigned int i = 0; i < nFun; i++){
fullVector<double> vi = fullVector<double> vi =
...@@ -57,5 +57,6 @@ interpolate(const MElement& element, ...@@ -57,5 +57,6 @@ interpolate(const MElement& element,
val.axpy(vi, 1); val.axpy(vi, 1);
} }
// Return Interpolated Value //
return val; return val;
} }
#include <vector>
#include "BasisScalar.h"
#include "Mapper.h"
#include "FunctionSpaceNode.h" #include "FunctionSpaceNode.h"
using namespace std;
FunctionSpaceNode::FunctionSpaceNode(const GroupOfElement& goe, FunctionSpaceNode::FunctionSpaceNode(const GroupOfElement& goe,
int order){ int order){
// Build 0Form Basis // // Build 0Form Basis //
...@@ -14,5 +20,36 @@ interpolate(const MElement& element, ...@@ -14,5 +20,36 @@ interpolate(const MElement& element,
const std::vector<double>& coef, const std::vector<double>& coef,
const fullVector<double>& xyz) const{ const fullVector<double>& xyz) const{
return 42; // Const Cast For MElement //
MElement& eelement =
const_cast<MElement&>(element);
// Get Basis Functions //
const vector<Polynomial>& fun = getBasis(element).getFunctions();
const unsigned int nFun = fun.size();
// Get Jacobian //
fullMatrix<double> invJac(3, 3);
eelement.getJacobian(xyz(0), xyz(1), xyz(2), invJac);
invJac.invertInPlace();
// Get First Vertex of Element (Physical Space) //
fullVector<double> origin(3);
origin(0) = eelement.getVertex(0)->x();
origin(1) = eelement.getVertex(0)->y();
origin(2) = eelement.getVertex(0)->z();
// Map XYZ to Reference Plane (UVW) //
fullVector<double> uvw = Mapper::invMap(xyz, origin, invJac);
// 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;
} }
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
#define _FUNCTIONSPACESCALAR_H_ #define _FUNCTIONSPACESCALAR_H_
#include "fullMatrix.h" #include "fullMatrix.h"
#include "BasisScalar.h"
#include "FunctionSpace.h" #include "FunctionSpace.h"
class FunctionSpaceScalar : public FunctionSpace{ class FunctionSpaceScalar : public FunctionSpace{
...@@ -12,6 +13,18 @@ class FunctionSpaceScalar : public FunctionSpace{ ...@@ -12,6 +13,18 @@ class FunctionSpaceScalar : public FunctionSpace{
interpolate(const MElement& element, interpolate(const MElement& element,
const std::vector<double>& coef, const std::vector<double>& coef,
const fullVector<double>& xyz) const = 0; const fullVector<double>& xyz) const = 0;
const BasisScalar& getBasis(const MElement& element) const;
}; };
//////////////////////
// Inline Functions //
//////////////////////
inline const BasisScalar& FunctionSpaceScalar::
getBasis(const MElement& element) const{
return static_cast<const BasisScalar&>(*basis);
}
#endif #endif
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
#define _FUNCTIONSPACEVECTOR_H_ #define _FUNCTIONSPACEVECTOR_H_
#include "fullMatrix.h" #include "fullMatrix.h"
#include "BasisVector.h"
#include "FunctionSpace.h" #include "FunctionSpace.h"
class FunctionSpaceVector : public FunctionSpace{ class FunctionSpaceVector : public FunctionSpace{
...@@ -12,6 +13,18 @@ class FunctionSpaceVector : public FunctionSpace{ ...@@ -12,6 +13,18 @@ class FunctionSpaceVector : public FunctionSpace{
interpolate(const MElement& element, interpolate(const MElement& element,
const std::vector<double>& coef, const std::vector<double>& coef,
const fullVector<double>& xyz) const = 0; const fullVector<double>& xyz) const = 0;
const BasisVector& getBasis(const MElement& element) const;
}; };
//////////////////////
// Inline Functions //
//////////////////////
inline const BasisVector& FunctionSpaceVector::
getBasis(const MElement& element) const{
return static_cast<const BasisVector&>(*basis);
}
#endif #endif
#include <cmath> #include <cmath>
#include <sstream> #include <sstream>
#include <stack> #include <stack>
#include "Polynomial.h" #include "Polynomial.h"
using namespace std; using namespace std;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment