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

Prepatation for Lagrange ReferenceSpace + correct consteness + cleaning some...

Prepatation for Lagrange ReferenceSpace + correct consteness + cleaning some tabs -- replaced by spaces
parent 6ba9a24a
Branches
Tags
No related merge requests found
......@@ -25,7 +25,7 @@ BasisHierarchical0From::~BasisHierarchical0From(void){
if(hasGrad){
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
delete grad[i][j];
delete grad[i][j];
delete[] grad[i];
}
......
#include "Exception.h"
#include "BasisLagrange.h"
using namespace std;
BasisLagrange::BasisLagrange(void){
scalar = true;
}
......@@ -10,30 +12,40 @@ BasisLagrange::~BasisLagrange(void){
unsigned int BasisLagrange::
getNOrientation(void) const{
throw Exception("BasisLagrange::Not Implemented");
return refSpace->getNPermutation();
}
unsigned int BasisLagrange::
getOrientation(const MElement& element) const{
throw Exception("BasisLagrange::Not Implemented");
return refSpace->getPermutation(element);
}
fullMatrix<double>* BasisLagrange::
getFunctions(const MElement& element,
double u, double v, double w) const{
throw Exception("BasisLagrange::Not Implemented");
}
double u, double v, double w) const{
fullMatrix<double>* BasisLagrange::getFunctions(unsigned int orientation,
double u, double v, double w) const{
// Fill Matrix //
fullMatrix<double> tmp;
fullMatrix<double> point(1, 3);
point(0, 0) = u;
point(0, 1) = v;
point(0, 2) = w;
throw Exception("BasisLagrange::Not Implemented");
/*
// Alloc Memory //
fullMatrix<double> tmp(nFunction, 1);
fullMatrix<double>* values = new fullMatrix<double>(nFunction, 1);
lBasis->f(point, tmp);
// Transpose 'tmp': otherwise not coherent with df !!
fullMatrix<double> values = tmp.transpose();
// Get Inorder & Return //
return inorder(refSpace->getPermutation(element), values);
}
fullMatrix<double>* BasisLagrange::
getFunctions(unsigned int orientation,
double u, double v, double w) const{
// Fill Matrix //
fullMatrix<double> tmp;
fullMatrix<double> point(1, 3);
point(0, 0) = u;
point(0, 1) = v;
......@@ -41,18 +53,19 @@ fullMatrix<double>* BasisLagrange::getFunctions(unsigned int orientation,
lBasis->f(point, tmp);
*values = tmp.transpose(); // Otherwise not coherent with df !!
// Transpose 'tmp': otherwise not coherent with df !!
fullMatrix<double> values = tmp.transpose();
// Return //
return values;
*/
// Get Inorder & Return //
return inorder(orientation, values);
}
void BasisLagrange::preEvaluateFunctions(const fullMatrix<double>& point) const{
throw Exception("BasisLagrange::Not Implemented");
}
void BasisLagrange::preEvaluateDerivatives(const fullMatrix<double>& point) const{
void BasisLagrange::
preEvaluateDerivatives(const fullMatrix<double>& point) const{
throw Exception("BasisLagrange::Not Implemented");
}
......@@ -76,14 +89,14 @@ BasisLagrange::getPreEvaluatedDerivatives(unsigned int orientation) const{
throw Exception("BasisLagrange::Not Implemented");
}
std::vector<double> BasisLagrange::
vector<double> BasisLagrange::
project(const MElement& element,
const std::vector<double>& coef,
const FunctionSpaceScalar& fSpace){
const vector<double>& coef,
const FunctionSpaceScalar& fSpace){
// Init New Coefs //
const unsigned int size = lPoint->size1();
std::vector<double> newCoef(size);
vector<double> newCoef(size);
// Interpolation at Lagrange Points //
for(unsigned int i = 0; i < size; i++){
......@@ -93,22 +106,22 @@ project(const MElement& element,
uvw(2) = (*lPoint)(i, 2);
newCoef[i] = fSpace.interpolateInRefSpace(element,
coef,
uvw);
coef,
uvw);
}
// Return ;
return newCoef;
}
std::vector<fullVector<double> > BasisLagrange::
vector<fullVector<double> > BasisLagrange::
project(const MElement& element,
const std::vector<double>& coef,
const FunctionSpaceVector& fSpace){
const vector<double>& coef,
const FunctionSpaceVector& fSpace){
// Init New Coefs //
const unsigned int size = lPoint->size1();
std::vector<fullVector<double> > newCoef(size);
vector<fullVector<double> > newCoef(size);
// Interpolation at Lagrange Points //
for(unsigned int i = 0; i < size; i++){
......@@ -118,10 +131,32 @@ project(const MElement& element,
uvw(2) = (*lPoint)(i, 2);
newCoef[i] = fSpace.interpolateInRefSpace(element,
coef,
uvw);
coef,
uvw);
}
// Return ;
return newCoef;
}
fullMatrix<double>* BasisLagrange::inorder(unsigned int orientation,
fullMatrix<double>& mat) const{
// Matrix Size //
const int nRow = mat.size1();
const int nCol = mat.size2();
// Order //
const vector<unsigned int>* order =
refSpace->getAllLagrangeNode()[orientation];
// Alloc new matrix //
fullMatrix<double>* ret = new fullMatrix<double>(nRow, nCol);
// Populate //
for(int i = 0; i < nRow; i++)
for(int j = 0; j < nCol; j++)
(*ret)(i, j) = mat((*order)[i], j);
// Return //
return ret;
}
......@@ -6,6 +6,7 @@
#include "FunctionSpaceVector.h"
#include "fullMatrix.h"
#include "polynomialBasis.h"
#include "ReferenceSpaceLagrange.h"
/**
@interface BasisLagrange
......@@ -25,8 +26,9 @@
class BasisLagrange: public BasisLocal{
protected:
polynomialBasis* lBasis; // Lagrange Basis
fullMatrix<double>* lPoint; // Lagrange Points
polynomialBasis* lBasis; // Lagrange Basis
fullMatrix<double>* lPoint; // Lagrange Points
ReferenceSpaceLagrange* refSpace; // Reference Space
public:
virtual ~BasisLagrange(void);
......@@ -35,10 +37,10 @@ class BasisLagrange: public BasisLocal{
virtual unsigned int getOrientation(const MElement& element) const;
virtual fullMatrix<double>* getFunctions(const MElement& element,
double u, double v, double w) const;
double u, double v, double w) const;
virtual fullMatrix<double>* getFunctions(unsigned int orientation,
double u, double v, double w) const;
double u, double v, double w) const;
virtual void preEvaluateFunctions(const fullMatrix<double>& point) const;
virtual void preEvaluateDerivatives(const fullMatrix<double>& point) const;
......@@ -60,16 +62,19 @@ class BasisLagrange: public BasisLocal{
std::vector<double>
project(const MElement& element,
const std::vector<double>& coef,
const FunctionSpaceScalar& fSpace);
const std::vector<double>& coef,
const FunctionSpaceScalar& fSpace);
std::vector<fullVector<double> >
project(const MElement& element,
const std::vector<double>& coef,
const FunctionSpaceVector& fSpace);
const std::vector<double>& coef,
const FunctionSpaceVector& fSpace);
protected:
BasisLagrange(void);
fullMatrix<double>* inorder(unsigned int orientation,
fullMatrix<double>& mat) const;
};
......
......@@ -13,6 +13,9 @@ set(SRC
QuadReferenceSpace.cpp
TetReferenceSpace.cpp
ReferenceSpaceLagrange.cpp
TriLagrangeReferenceSpace.cpp
Basis.cpp
BasisLocal.cpp
BasisGenerator.cpp
......
......@@ -158,17 +158,17 @@ void ReferenceSpace::populate(node* pTreeRoot){
pTreeRoot->next[i].last[depth] = nextLast;
for(unsigned int j = 0; j < depth; j++)
pTreeRoot->next[i].last[j] = pTreeRoot->last[j];
pTreeRoot->next[i].last[j] = pTreeRoot->last[j];
// Possibilities of child node
pTreeRoot->next[i].number = nextNumber;
pTreeRoot->next[i].possible = new unsigned int[nextNumber];
for(unsigned int j = 0; j < nextNumber; j++){
if(pTreeRoot->possible[j] == nextLast)
offset = 1;
if(pTreeRoot->possible[j] == nextLast)
offset = 1;
pTreeRoot->next[i].possible[j] = pTreeRoot->possible[j + offset];
pTreeRoot->next[i].possible[j] = pTreeRoot->possible[j + offset];
}
// Populate each child node (until a leaf is found)
......@@ -199,8 +199,8 @@ void ReferenceSpace::getEdge(void){
for(unsigned int i = 0; i < nEdge; i++)
(*tmp)[i] = inOrder(p,
refEdge[i][0],
refEdge[i][1]);
refEdge[i][0],
refEdge[i][1]);
(*edge)[p] = tmp;
}
}
......@@ -214,9 +214,9 @@ void ReferenceSpace::getTriFace(void){
for(unsigned int i = 0; i < nFace; i++)
(*tmp)[i] = inOrder(p,
refFace[i][0],
refFace[i][1],
refFace[i][2]);
refFace[i][0],
refFace[i][1],
refFace[i][2]);
(*face)[p] = tmp;
}
}
......@@ -230,9 +230,9 @@ void ReferenceSpace::getQuaFace(void){
for(unsigned int i = 0; i < nFace; i++)
(*tmp)[i] = inOrder(p,
refFace[i][0],
refFace[i][1],
refFace[i][2],
refFace[i][0],
refFace[i][1],
refFace[i][2],
refFace[i][3]);
(*face)[p] = tmp;
}
......@@ -240,8 +240,8 @@ void ReferenceSpace::getQuaFace(void){
const vector<unsigned int>* ReferenceSpace::
inOrder(unsigned int permutation,
unsigned int a,
unsigned int b){
unsigned int a,
unsigned int b){
unsigned int v;
unsigned int k = 0;
......@@ -262,9 +262,9 @@ inOrder(unsigned int permutation,
const std::vector<unsigned int>* ReferenceSpace::
inOrder(unsigned int permutation,
unsigned int a,
unsigned int b,
unsigned int c){
unsigned int a,
unsigned int b,
unsigned int c){
unsigned int v;
unsigned int k = 0;
......@@ -285,9 +285,9 @@ inOrder(unsigned int permutation,
const std::vector<unsigned int>* ReferenceSpace::
inOrder(unsigned int permutation,
unsigned int a,
unsigned int b,
unsigned int c,
unsigned int a,
unsigned int b,
unsigned int c,
unsigned int d){
unsigned int opposite;
......@@ -354,13 +354,13 @@ unsigned int ReferenceSpace::getPermutation(const MElement& elem) const{
catch(...){
throw Exception("Cannot Find Reference Space for Element %d",
element.getNum());
element.getNum());
}
}
unsigned int ReferenceSpace::treeLookup(const node* root,
vector<pair<unsigned int, MVertex*> >&
sortedArray){
vector<pair<unsigned int, MVertex*> >&
sortedArray){
// Temp Data //
unsigned int choice;
unsigned int i;
......@@ -414,8 +414,8 @@ string ReferenceSpace::toString(void) const{
for(unsigned int j = 0; j < nEdge; j++)
stream << " -- ["
<< edge->at(i)->at(j)->at(0) << ", "
<< edge->at(i)->at(j)->at(1) << "]" << endl;
<< edge->at(i)->at(j)->at(0) << ", "
<< edge->at(i)->at(j)->at(1) << "]" << endl;
}
stream << "Faces Permutations:" << endl;
......@@ -425,9 +425,9 @@ string ReferenceSpace::toString(void) const{
for(unsigned int j = 0; j < nFace; j++)
stream << " -- ["
<< face->at(i)->at(j)->at(0) << ", "
<< face->at(i)->at(j)->at(1) << ", "
<< face->at(i)->at(j)->at(2) << "]" << endl;
<< face->at(i)->at(j)->at(0) << ", "
<< face->at(i)->at(j)->at(1) << ", "
<< face->at(i)->at(j)->at(2) << "]" << endl;
}
return stream.str();
......@@ -449,3 +449,16 @@ string ReferenceSpace::toString(const node* node) const{
return stream.str();
}
string ReferenceSpace::toLatex(void) const{
stringstream stream;
stream << "\\documentclass{article}" << endl << endl
<< "\\begin{document}" << endl
<< "\texttt{toLatex} not implemented" << endl
<< "\\end{document}" << endl;
return stream.str();
}
......@@ -67,15 +67,15 @@ class ReferenceSpace{
unsigned int getPermutation(const MElement& element) const;
std::vector<const std::vector<const std::vector<unsigned int>*>*>&
const std::vector<const std::vector<const std::vector<unsigned int>*>*>&
getAllEdge(void) const;
std::vector<const std::vector<const std::vector<unsigned int>*>*>&
const std::vector<const std::vector<const std::vector<unsigned int>*>*>&
getAllFace(void) const ;
std::string toString(void) const;
virtual std::string toString(void) const;
virtual std::string toLatex(void) const = 0;
virtual std::string toLatex(void) const;
protected:
ReferenceSpace(void);
......@@ -93,26 +93,26 @@ class ReferenceSpace{
void getQuaFace(void);
const std::vector<unsigned int>* inOrder(unsigned int permutation,
unsigned int a,
unsigned int b);
unsigned int a,
unsigned int b);
const std::vector<unsigned int>* inOrder(unsigned int permutation,
unsigned int a,
unsigned int b,
unsigned int c);
unsigned int a,
unsigned int b,
unsigned int c);
const std::vector<unsigned int>* inOrder(unsigned int permutation,
unsigned int a,
unsigned int b,
unsigned int c,
unsigned int a,
unsigned int b,
unsigned int c,
unsigned int d);
static bool sortPredicate(const std::pair<unsigned int, MVertex*>& a,
const std::pair<unsigned int, MVertex*>& b);
const std::pair<unsigned int, MVertex*>& b);
static unsigned int treeLookup(const node* root,
std::vector<std::pair<unsigned int, MVertex*> >&
sortedArray);
std::vector<std::pair<unsigned int, MVertex*> >&
sortedArray);
std::string toString(const node* node) const;
};
......@@ -188,13 +188,13 @@ ReferenceSpace::getNPermutation(void) const{
}
inline
std::vector<const std::vector<const std::vector<unsigned int>*>*>&
const std::vector<const std::vector<const std::vector<unsigned int>*>*>&
ReferenceSpace::getAllEdge(void) const{
return *edge;
}
inline
std::vector<const std::vector<const std::vector<unsigned int>*>*>&
const std::vector<const std::vector<const std::vector<unsigned int>*>*>&
ReferenceSpace::getAllFace(void) const{
return *face;
}
......@@ -202,7 +202,7 @@ ReferenceSpace::getAllFace(void) const{
inline
bool
ReferenceSpace::sortPredicate(const std::pair<unsigned int, MVertex*>& a,
const std::pair<unsigned int, MVertex*>& b){
const std::pair<unsigned int, MVertex*>& b){
return a.second->getNum() < b.second->getNum();
}
......
#include <sstream>
#include "ReferenceSpaceLagrange.h"
using namespace std;
ReferenceSpaceLagrange::ReferenceSpaceLagrange(void){
// Init to NULL //
node = NULL;
}
ReferenceSpaceLagrange::~ReferenceSpaceLagrange(void){
// If 'node' allocated //
if(node){
for(unsigned int p = 0; p < nPerm; p++)
delete (*node)[p];
delete node;
}
}
void ReferenceSpaceLagrange::getLagrangeNode(void){
// Alloc //
vector<unsigned int>* tmp;
node = new vector<const vector<unsigned int>*>(nPerm);
// Populate //
for(unsigned int p = 0; p < nPerm; p++){
// Alloc Temp //
tmp = new vector<unsigned int>(nNode);
// Vertex based node
for(unsigned int i = 0; i < nVertex; i++)
(*tmp)[i] = i;
// Edge based node
for(unsigned int e = 0; e < nEdge; e++)
edgeSeq(*tmp,
nVertex + nNodePerEdge * e,
nVertex + nNodePerEdge * e,
nVertex + nNodePerEdge * (e + 1),
refEdge[e],
*(*(*edge)[p])[e]);
// Face based node
for(unsigned int f = 0; f < nFace; f++)
faceSeq(*tmp,
nVertex + nNodePerEdge * nEdge + nNodePerFace * f,
nVertex + nNodePerEdge * nEdge + nNodePerFace * f,
nVertex + nNodePerEdge * nEdge + nNodePerFace * (f + 1),
refFace[f],
*(*(*face)[p])[f],
nNodePerEdge);
// Insert in node
(*node)[p] = tmp;
}
}
void ReferenceSpaceLagrange::edgeSeq(vector<unsigned int>& vec,
unsigned int startIdx,
unsigned int startVal,
unsigned int stopVal,
unsigned int* refEdge,
const vector<unsigned int>& edge){
// Is reverted ? //
const bool isRevert = (edge[0] != refEdge[0]);
// Index //
unsigned int val = startVal;
unsigned int idx;
// Depending if edge is reverted ...
if(isRevert)
idx = startIdx + stopVal - startVal - 1;
else
idx = startIdx;
// Populate //
for(; val < stopVal; val++){
vec[idx] = val;
if(isRevert)
idx--;
else
idx++;
}
}
void ReferenceSpaceLagrange::faceSeq(vector<unsigned int>& vec,
unsigned int startIdx,
unsigned int startVal,
unsigned int stopVal,
unsigned int* refFace,
const vector<unsigned int>& face,
unsigned int nNodePerEdge){
// Index //
unsigned int val = startVal;
unsigned int idx = startIdx;
// Populate //
for(; val < stopVal; val++){
vec[idx] = val;
idx++;
}
}
string ReferenceSpaceLagrange::toString(void) const{
const unsigned int nNodeMinus = nNode - 1;
stringstream stream;
stream << ReferenceSpace::toString();
stream << "Lagrange Nodes Permutations:" << endl;
for(unsigned int i = 0; i < nPerm; i++){
stream << " * RefSpace #" << i + 1 << ":" << endl
<< " -- [";
for(unsigned int j = 0; j < nNodeMinus; j++)
stream << node->at(i)->at(j) << ", ";
stream << node->at(i)->at(nNodeMinus) << "]" << endl;
}
return stream.str();
}
#ifndef _REFERENCESPACELAGRANGE_H_
#define _REFERENCESPACELAGRANGE_H_
#include "ReferenceSpace.h"
/**
@interface ReferenceSpaceLagrange
@brief Base interface for all Lagrange ReferenceSpace%s
This class is a particular ReferenceSpace for
@em Lagrange Basis.@n
This class adds the notion of Lagrange points and
their permutations depending on the orientation
of a MElement.@n
Lagrange points are ordered in a way such that
Lagrange function must be taken in that order in order
to have an oriented basis.@n
Because the number of Lagrange Points depends on
the the order of the Lagrange Basis,
a ReferenceSpaceLagrange shall depend on the order.
*/
class ReferenceSpaceLagrange: public ReferenceSpace{
protected:
// Lagrange Node Permutation //
unsigned int nNode;
unsigned int nNodePerEdge;
unsigned int nNodePerFace;
unsigned int nNodePerCell;
std::vector<const std::vector<unsigned int>*>* node;
public:
virtual ~ReferenceSpaceLagrange(void);
const std::vector<const std::vector<unsigned int>*>&
getAllLagrangeNode(void) const;
virtual std::string toString(void) const;
protected:
ReferenceSpaceLagrange(void);
void getLagrangeNode(void);
static void edgeSeq(std::vector<unsigned int>& vec,
unsigned int startIdx,
unsigned int startVal,
unsigned int stopVal,
unsigned int* refEdge,
const std::vector<unsigned int>& edge);
static void faceSeq(std::vector<unsigned int>& vec,
unsigned int startIdx,
unsigned int startVal,
unsigned int stopVal,
unsigned int* refFace,
const std::vector<unsigned int>& face,
unsigned int nNodePerEdge);
};
/////////////////////
// Inline Function //
/////////////////////
inline
const std::vector<const std::vector<unsigned int>*>&
ReferenceSpaceLagrange::getAllLagrangeNode(void) const{
return *node;
}
#endif
#include "Exception.h"
#include "TriLagrangeBasis.h"
#include "pointsGenerators.h"
#include "TriLagrangeReferenceSpace.h"
#include <iostream>
TriLagrangeBasis::TriLagrangeBasis(unsigned int order){
// If order 0 (Nedelec): use order 1
......@@ -25,11 +28,16 @@ TriLagrangeBasis::TriLagrangeBasis(unsigned int order){
// Init Lagrange Point //
lPoint = new fullMatrix<double>
(gmshGeneratePointsTriangle(order, false));
// ReferenceSpace //
refSpace = new TriLagrangeReferenceSpace(order);
std::cout << refSpace->toString() << std::endl;
}
TriLagrangeBasis::~TriLagrangeBasis(void){
delete lBasis;
delete lPoint;
delete refSpace;
}
unsigned int TriLagrangeBasis::getTag(unsigned int order){
......
......@@ -50,30 +50,30 @@ string TriReferenceSpace::toLatex(void) const{
stream << "\\documentclass{article}" << endl << endl
<< "\\usepackage{longtable}" << endl
<< "\\usepackage{tikz}" << endl
<< "\\usetikzlibrary{arrows}" << endl << endl
<< "\\usepackage{longtable}" << endl
<< "\\usepackage{tikz}" << endl
<< "\\usetikzlibrary{arrows}" << endl << endl
<< "\\begin{document}" << endl
<< "\\tikzstyle{vertex} = [circle, fill = black!25]" << endl
<< "\\tikzstyle{line} = [draw, thick, black, -latex']" << endl << endl
<< "\\begin{document}" << endl
<< "\\tikzstyle{vertex} = [circle, fill = black!25]" << endl
<< "\\tikzstyle{line} = [draw, thick, black, -latex']" << endl << endl
<< "\\begin{longtable}{ccc}" << endl << endl;
<< "\\begin{longtable}{ccc}" << endl << endl;
for(unsigned int p = 0; p < nPerm; p++){
stream << "\\begin{tikzpicture}" << endl
<< "\\node[vertex] (n0) at(0, 0) {$" << perm[p][0] << "$};" << endl
<< "\\node[vertex] (n1) at(3, 0) {$" << perm[p][1] << "$};" << endl
<< "\\node[vertex] (n2) at(0, 3) {$" << perm[p][2] << "$};" << endl
<< "\\node[vertex] (n0) at(0, 0) {$" << perm[p][0] << "$};" << endl
<< "\\node[vertex] (n1) at(3, 0) {$" << perm[p][1] << "$};" << endl
<< "\\node[vertex] (n2) at(0, 3) {$" << perm[p][2] << "$};" << endl
<< endl;
for(unsigned int i = 0; i < 3; i++)
stream << "\\path[line]"
<< " (n" << (*(*(*edge)[p])[i])[0] << ")"
<< " -- "
<< " (n" << (*(*(*edge)[p])[i])[1] << ");"
<< endl;
<< " (n" << (*(*(*edge)[p])[i])[0] << ")"
<< " -- "
<< " (n" << (*(*(*edge)[p])[i])[1] << ");"
<< endl;
if((p + 1) % 3)
stream << "\\end{tikzpicture} & " << endl << endl;
......@@ -83,7 +83,7 @@ string TriReferenceSpace::toLatex(void) const{
}
stream << "\\end{longtable}" << endl
<< "\\end{document}" << endl;
<< "\\end{document}" << endl;
return stream.str();
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment