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

ReferenceSpace computes Jacobians: Seems OK -- TriEdgeBasis seems to be OK --... 

ReferenceSpace computes Jacobians: Seems OK -- TriEdgeBasis seems to be OK -- FunctionSpaceScalar and FunctionSpaceVector adapted with new ReferenceSpace::getJacobian
parent 9d60f13e
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 368 additions and 205 deletions
FunctionSpace/BasisGenerator.cpp
+ 17
− 15
View file @ f529232c
......@@ -147,6 +147,8 @@ BasisLocal* BasisGenerator::hexHierarchicalGen(unsigned int basisType,
switch(basisType){
//case 0: return new HexNodeBasis(order);
//case 1: return new HexEdgeBasis(order);
case 0: throw Exception("0-form not implemented on Hexs");
case 1: throw Exception("1-form not implemented on Hexs");
case 2: throw Exception("2-form not implemented on Hexs");
case 3: throw Exception("3-form not implemented on Hexs");
......
FunctionSpace/FunctionSpaceScalar.cpp
+ 3
− 32
View file @ f529232c
......@@ -11,38 +11,9 @@ FunctionSpaceScalar::~FunctionSpaceScalar(void){
}
double FunctionSpaceScalar::
interpolate(const MElement& element,
interpolateInABC(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const{
// Get ABC Space coordinate //
double abc[3];
(*basis)[0]->getReferenceSpace().mapFromXYZtoABC(element, xyz, abc);
// Get Basis Functions //
const unsigned int nFun = (*basis)[0]->getNFunction();
fullMatrix<double> fun(nFun, 1);
(*basis)[0]->getFunctions(fun, element, abc[0], abc[1], abc[2]);
// Interpolate (in Reference Place) //
double val = 0;
for(unsigned int i = 0; i < nFun; i++)
val += fun(i, 0) * coef[i];
// Return Interpolated Value //
return val;
}
double FunctionSpaceScalar::
interpolateInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const{
// Get ABC Space coordinate //
double abc[3];
(*basis)[0]->getReferenceSpace().mapFromUVWtoABC(element, uvw, abc);
double abc[3]) const{
// Get Basis Functions //
const unsigned int nFun = (*basis)[0]->getNFunction();
......
FunctionSpace/FunctionSpaceScalar.h
+ 38
− 0
View file @ f529232c
......@@ -30,6 +30,11 @@ class FunctionSpaceScalar : public FunctionSpace{
interpolateInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const;
private:
double interpolateInABC(const MElement& element,
const std::vector<double>& coef,
double abc[3]) const;
};
......@@ -84,4 +89,37 @@ class FunctionSpaceScalar : public FunctionSpace{
---> check
*/
//////////////////////
// Inline Functions //
//////////////////////
inline double FunctionSpaceScalar::
interpolate(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const{
// Get ABC Space coordinate //
double abc[3];
(*basis)[0]->getReferenceSpace().mapFromXYZtoABC(element,
xyz(0), xyz(1), xyz(2),
abc);
// Interpolate in ABC //
return interpolateInABC(element, coef, abc);
}
inline double FunctionSpaceScalar::
interpolateInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const{
// Get ABC Space coordinate //
double abc[3];
(*basis)[0]->getReferenceSpace().mapFromUVWtoABC(element,
uvw(0), uvw(1), uvw(2),
abc);
// Interpolate in ABC //
return interpolateInABC(element, coef, abc);
}
#endif
FunctionSpace/FunctionSpaceVector.cpp
+ 7
− 53
View file @ f529232c
......@@ -12,68 +12,22 @@ FunctionSpaceVector::~FunctionSpaceVector(void){
}
fullVector<double> FunctionSpaceVector::
interpolate(const MElement& element,
interpolateInABC(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const{
// Const Cast For MElement //
MElement& eelement =
const_cast<MElement&>(element);
// Get Reference coordinate //
double phys[3] = {xyz(0), xyz(1), xyz(2)};
double uvw[3];
eelement.xyz2uvw(phys, uvw);
// Get Jacobian //
fullMatrix<double> invJac(3, 3);
eelement.getJacobian(uvw[0], uvw[1], uvw[2], invJac);
invJac.invertInPlace();
// Get Basis Functions //
const unsigned int nFun = (*basis)[0]->getNFunction();
fullMatrix<double> fun(nFun, 3);
(*basis)[0]->getFunctions(fun, element, uvw[0], uvw[1], uvw[2]);
// Interpolate (in Reference Place) //
fullMatrix<double> val(1, 3);
val(0, 0) = 0;
val(0, 1) = 0;
val(0, 2) = 0;
for(unsigned int i = 0; i < nFun; i++){
val(0, 0) += fun(i, 0) * coef[i];
val(0, 1) += fun(i, 1) * coef[i];
val(0, 2) += fun(i, 2) * coef[i];
}
// Return Interpolated Value //
fullVector<double> map(3);
Mapper::hCurl(val, 0, 0, invJac, map);
return map;
}
fullVector<double> FunctionSpaceVector::
interpolateInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const{
// Const Cast For MElement //
MElement& eelement =
const_cast<MElement&>(element);
double abc[3]) const{
// Get Jacobian //
fullMatrix<double> invJac(3, 3);
eelement.getJacobian(uvw(0), uvw(1), uvw(2), invJac);
(*basis)[0]->getReferenceSpace().getJacobian(element,
abc[0], abc[1], abc[2],
invJac);
invJac.invertInPlace();
// Get Basis Functions //
const unsigned int nFun = (*basis)[0]->getNFunction();
fullMatrix<double> fun(nFun, 3);
(*basis)[0]->getFunctions(fun, element, uvw(0), uvw(1), uvw(2));
(*basis)[0]->getFunctions(fun, element, abc[0], abc[1], abc[2]);
// Interpolate (in Reference Place) //
fullMatrix<double> val(1, 3);
......
FunctionSpace/FunctionSpaceVector.h
+ 39
− 0
View file @ f529232c
......@@ -31,6 +31,12 @@ class FunctionSpaceVector : public FunctionSpace{
interpolateInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const;
private:
fullVector<double>
interpolateInABC(const MElement& element,
const std::vector<double>& coef,
double abc[3]) const;
};
......@@ -85,4 +91,37 @@ class FunctionSpaceVector : public FunctionSpace{
---> check
*/
//////////////////////
// Inline Functions //
//////////////////////
inline fullVector<double> FunctionSpaceVector::
interpolate(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& xyz) const{
// Get ABC Space coordinate //
double abc[3];
(*basis)[0]->getReferenceSpace().mapFromXYZtoABC(element,
xyz(0), xyz(1), xyz(2),
abc);
// Interpolate in ABC //
return interpolateInABC(element, coef, abc);
}
inline fullVector<double> FunctionSpaceVector::
interpolateInRefSpace(const MElement& element,
const std::vector<double>& coef,
const fullVector<double>& uvw) const{
// Get ABC Space coordinate //
double abc[3];
(*basis)[0]->getReferenceSpace().mapFromUVWtoABC(element,
uvw(0), uvw(1), uvw(2),
abc);
// Interpolate in ABC //
return interpolateInABC(element, coef, abc);
}
#endif
FunctionSpace/ReferenceSpace.cpp
+ 190
− 27
View file @ f529232c
......@@ -2,6 +2,7 @@
#include <sstream>
#include "Exception.h"
#include "Numeric.h"
#include "ReferenceSpace.h"
using namespace std;
......@@ -419,69 +420,231 @@ size_t ReferenceSpace::getPermutationIdx(const MElement& element) const{
}
}
void ReferenceSpace::mapFromABCtoUVW(const MElement& element,
double a, double b, double c,
double uvw[3]) const{
// Get Index Permutation
const size_t permutationIdx = getPermutationIdx(element);
// UVW node coordinate
double** uvwNode = new double*[nVertex];
for(size_t i = 0; i < nVertex; i++)
uvwNode[i] = new double[3];
for(size_t i = 0; i < nVertex; i++)
element.getNode(i,
uvwNode[UVWtoABCIndex[permutationIdx][i]][0],
uvwNode[UVWtoABCIndex[permutationIdx][i]][1],
uvwNode[UVWtoABCIndex[permutationIdx][0]][2]);
// ABC (order 1) grad shape functions
double* phiABC = new double[nVertex];
element.getShapeFunctions(a, b, c, phiABC, 1);
// Map From UVW to UVW //
uvw[0] = 0;
for(size_t i = 0; i < nVertex; i++)
uvw[0] += uvwNode[i][0] * phiABC[i];
uvw[1] = 0;
for(size_t i = 0; i < nVertex; i++)
uvw[1] += uvwNode[i][1] * phiABC[i];
uvw[2] = 0;
for(size_t i = 0; i < nVertex; i++)
uvw[2] += uvwNode[i][2] * phiABC[i];
// Free //
delete[] phiABC;
for(size_t i = 0; i < nVertex; i++)
delete[] uvwNode[i];
delete[] uvwNode;
}
void ReferenceSpace::mapFromXYZtoABC(const MElement& element,
const fullVector<double>& xyz,
double x, double y, double z,
double abc[3]) const{
// Get UVW coordinate //
double phys[3] = {xyz(0), xyz(1), xyz(2)};
double xyz[3] = {x, y, z};
double uvw[3];
element.xyz2uvw(phys, uvw);
element.xyz2uvw(xyz, uvw);
// Get ABC coordinate //
fullVector<double> uuvw(3);
uuvw(0) = uvw[0];
uuvw(1) = uvw[1];
uuvw(2) = uvw[2];
mapFromUVWtoABC(element, uuvw, abc);
mapFromUVWtoABC(element, uvw[0], uvw[1], uvw[2], abc);
}
void ReferenceSpace::mapFromUVWtoABC(const MElement& element,
const fullVector<double>& uvw,
double u, double v, double w,
double abc[3]) const{
// Compute coordinate in ABC Space //
// ABC node coordinate
double** abcMat = new double*[nVertex];
double** abcNode = new double*[nVertex];
for(size_t i = 0; i < nVertex; i++)
abcMat[i] = new double[3];
abcNode[i] = new double[3];
for(size_t i = 0; i < nVertex; i++)
element.getNode(i, abcMat[i][0], abcMat[i][1], abcMat[i][2]);
element.getNode(i, abcNode[i][0], abcNode[i][1], abcNode[i][2]);
// UVW (order 1) shape functions
double* phiUVW = new double[nVertex];
element.getShapeFunctions(uvw(0), uvw(1), uvw(2), phiUVW, 1);
element.getShapeFunctions(u, v, w, phiUVW, 1);
// Element Permutation Index //
// Element Permutation Index
const size_t permutationIdx = getPermutationIdx(element);
// Map From UVW to ABC //
// Map From UVW to ABC
abc[0] = 0;
for(size_t i = 0; i < nVertex; i++)
abc[0] += abcMat[i][0] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
abc[0] += abcNode[i][0] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
abc[1] = 0;
for(size_t i = 0; i < nVertex; i++)
abc[1] += abcMat[i][1] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
abc[1] += abcNode[i][1] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
abc[2] = 0;
for(size_t i = 0; i < nVertex; i++)
abc[2] += abcMat[i][2] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
abc[2] += abcNode[i][2] * phiUVW[ABCtoUVWIndex[permutationIdx][i]];
// Free //
// Free
delete[] phiUVW;
for(size_t i = 0; i < nVertex; i++)
delete[] abcMat[i];
delete[] abcNode[i];
delete[] abcMat;
delete[] abcNode;
}
void ReferenceSpace::getJacobian(const MElement& element,
const fullVector<double>& xyz,
double ReferenceSpace::getJacobian(const MElement& element,
double a, double b, double c,
fullMatrix<double>& jac) const{
// ABC to UVW Jacobian //
// Get Index Permutation
const size_t permutationIdx = getPermutationIdx(element);
// UVW node coordinate
double** uvwNode = new double*[nVertex];
for(size_t i = 0; i < nVertex; i++)
uvwNode[i] = new double[3];
for(size_t i = 0; i < nVertex; i++)
element.getNode(i,
uvwNode[UVWtoABCIndex[permutationIdx][i]][0],
uvwNode[UVWtoABCIndex[permutationIdx][i]][1],
uvwNode[UVWtoABCIndex[permutationIdx][i]][2]);
// ABC (order 1) grad shape functions
double phiABC[1256][3]; // Cannot be dynamicaly allocated since
// GMSH wants a double[][3] for phiABC !!!
element.getGradShapeFunctions(a, b, c, phiABC, 1);
// Jacobian
fullMatrix<double> jacABCtoUVW(3, 3);
jacABCtoUVW.setAll(0);
for(size_t i = 0; i < nVertex; i++){
jacABCtoUVW(0, 0) += uvwNode[i][0] * phiABC[i][0];
jacABCtoUVW(0, 1) += uvwNode[i][1] * phiABC[i][0];
jacABCtoUVW(0, 2) += uvwNode[i][2] * phiABC[i][0];
}
for(size_t i = 0; i < nVertex; i++){
jacABCtoUVW(1, 0) += uvwNode[i][0] * phiABC[i][1];
jacABCtoUVW(1, 1) += uvwNode[i][1] * phiABC[i][1];
jacABCtoUVW(1, 2) += uvwNode[i][2] * phiABC[i][1];
}
for(size_t i = 0; i < nVertex; i++){
jacABCtoUVW(2, 0) += uvwNode[i][0] * phiABC[i][2];
jacABCtoUVW(2, 1) += uvwNode[i][1] * phiABC[i][2];
jacABCtoUVW(2, 2) += uvwNode[i][2] * phiABC[i][2];
}
// Regularize Jacobian
regularize(element.getDim(), jacABCtoUVW);
// Free
for(size_t i = 0; i < nVertex; i++)
delete[] uvwNode[i];
delete[] uvwNode;
// Map ABC coordinate to UVW point //
double uvw[3];
mapFromABCtoUVW(element, a, b, c, uvw);
// UVW to XYZ Jacobian + Determinant //
// NB: Volume is not modified when we go from ABC to UVW //
// --> Determinant unchanged //
fullMatrix<double> jacUVWtoXYZ(3, 3);
double det = element.getJacobian(uvw[0], uvw[1], uvw[2], jacUVWtoXYZ);
// Product of the two Jacobians & Return //
jac.gemm(jacABCtoUVW, jacUVWtoXYZ, 1, 0);
return det;
}
void ReferenceSpace::regularize(size_t dim, fullMatrix<double>& jac){
double a[3];
double b[3];
double c[3];
switch(dim){
case 0:
jac(0, 0) = jac(1, 1) = jac(2, 2) = 1.0;
jac(0, 1) = jac(1, 0) = jac(2, 0) = 0.0;
jac(0, 2) = jac(1, 2) = jac(2, 1) = 0.0;
break;
case 1:
a[0] = jac(0, 0);
a[1] = jac(0, 1);
a[2] = jac(0, 2);
if((fabs(a[0]) >= fabs(a[1]) && fabs(a[0]) >= fabs(a[2])) ||
(fabs(a[1]) >= fabs(a[0]) && fabs(a[1]) >= fabs(a[2]))){
b[0] = a[1];
b[1] = -a[0];
b[2] = 0.;
}
else{
b[0] = 0.;
b[1] = a[2];
b[2] = -a[1];
}
norme(b);
prodve(a, b, c);
norme(c);
jac(1, 0) = b[0]; jac(1, 1) = b[1]; jac(1, 2) = b[2];
jac(2, 0) = c[0]; jac(2, 1) = c[1]; jac(2, 2) = c[2];
break;
case 2:
a[0] = jac(0, 0);
a[1] = jac(0, 1);
a[2] = jac(0, 2);
b[0] = jac(1, 0);
b[1] = jac(1, 1);
b[2] = jac(1, 2);
prodve(a, b, c);
norme(c);
jac(2, 0) = c[0]; jac(2, 1) = c[1]; jac(2, 2) = c[2];
break;
case 3:
break;
}
}
string ReferenceSpace::toString(void) const{
......
FunctionSpace/ReferenceSpace.h
+ 11
− 5
View file @ f529232c
......@@ -64,16 +64,20 @@ class ReferenceSpace{
const std::vector<size_t>&
getNodeIndexFromABCtoUVW(const MElement& element) const;
void mapFromABCtoUVW(const MElement& element,
double a, double b, double c,
double uvw[3]) const;
void mapFromUVWtoABC(const MElement& element,
const fullVector<double>& xyz,
double u, double v, double w,
double abc[3]) const;
void mapFromXYZtoABC(const MElement& element,
const fullVector<double>& xyz,
double x, double y, double z,
double abc[3]) const;
void getJacobian(const MElement& element,
const fullVector<double>& xyz,
double getJacobian(const MElement& element,
double a, double b, double c,
fullMatrix<double>& jac) const;
virtual std::string toString(void) const;
......@@ -123,6 +127,8 @@ class ReferenceSpace{
static bool sortPredicate(const std::pair<size_t, size_t>& a,
const std::pair<size_t, size_t>& b);
static void regularize(size_t dim, fullMatrix<double>& jac);
};
......
FunctionSpace/TriEdgeBasis.cpp
+ 47
− 52
View file @ f529232c
......@@ -4,17 +4,16 @@
using namespace std;
TriEdgeBasis::TriEdgeBasis(unsigned int order){
/*
TriEdgeBasis::TriEdgeBasis(size_t order){
// Reference Space //
refSpace = new TriReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
const vector<vector<vector<size_t> > >&
edgeIdx = refSpace->getEdgeNodeIndex();
const vector<const vector<const vector<unsigned int>*>*>&
faceV = refSpace->getAllFace();
const vector<vector<vector<size_t> > >&
faceIdx = refSpace->getFaceNodeIndex();
// Set Basis Type //
this->order = order;
......@@ -54,28 +53,27 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
// Basis //
basis = new vector<Polynomial>**[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
for(size_t s = 0; s < nRefSpace; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based //
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = 0;
for(size_t s = 0; s < nRefSpace; s++){
size_t i = 0;
for(int e = 0; e < 3; e++){
for(int l = 0; l < orderPlus; l++){
// Nedelec
if(l == 0){
vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].gradient();
vector<Polynomial> tmp1 = lagrange[edgeIdx[s][e][1]].gradient();
vector<Polynomial> tmp2 = lagrange[edgeIdx[s][e][0]].gradient();
tmp1[0].mul(lagrange[(*(*edgeV[s])[e])[0]]);
tmp1[1].mul(lagrange[(*(*edgeV[s])[e])[0]]);
tmp1[2].mul(lagrange[(*(*edgeV[s])[e])[0]]);
tmp1[0].mul(lagrange[edgeIdx[s][e][0]]);
tmp1[1].mul(lagrange[edgeIdx[s][e][0]]);
tmp1[2].mul(lagrange[edgeIdx[s][e][0]]);
tmp2[0].mul(lagrange[(*(*edgeV[s])[e])[1]]);
tmp2[1].mul(lagrange[(*(*edgeV[s])[e])[1]]);
tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
tmp2[0].mul(lagrange[edgeIdx[s][e][1]]);
tmp2[1].mul(lagrange[edgeIdx[s][e][1]]);
tmp2[2].mul(lagrange[edgeIdx[s][e][1]]);
tmp2[0].sub(tmp1[0]);
tmp2[1].sub(tmp1[1]);
......@@ -88,11 +86,11 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
else{
basis[s][i] =
new vector<Polynomial>
((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
lagrange[(*(*edgeV[s])[e])[1]]
((intLegendre[l].compose(lagrange[edgeIdx[s][e][0]] -
lagrange[edgeIdx[s][e][1]]
,
lagrange[(*(*edgeV[s])[e])[1]] +
lagrange[(*(*edgeV[s])[e])[0]])).gradient());
lagrange[edgeIdx[s][e][1]] +
lagrange[edgeIdx[s][e][0]])).gradient());
}
i++;
}
......@@ -101,7 +99,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
// Face Based //
// NB: We use (*(*faceV[s])[f])[]
// NB: We use (*(*faceIdx[s])[f])[]
// where f = 0, because triangles
// have only ONE face: the face '0'
......@@ -113,40 +111,40 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
Polynomial* p = new Polynomial[nRefSpace];
vector<Polynomial>** subGrad = new vector<Polynomial>*[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
for(size_t s = 0; s < nRefSpace; s++)
u[s] = new Polynomial[orderPlus];
for(unsigned int s = 0; s < nRefSpace; s++)
for(size_t s = 0; s < nRefSpace; s++)
v[s] = new Polynomial[orderPlus];
// Preliminaries
for(unsigned int s = 0; s < nRefSpace; s++){
p[s] = lagrange[(*(*faceV[s])[0])[2]] * 2 - Polynomial(1, 0, 0, 0);
for(size_t s = 0; s < nRefSpace; s++){
p[s] = lagrange[faceIdx[s][0][2]] * 2 - Polynomial(1, 0, 0, 0);
// Polynomial u(x) & v(x)
for(int l = 0; l < orderPlus; l++){
u[s][l] = intLegendre[l].compose(lagrange[(*(*faceV[s])[0])[1]] -
lagrange[(*(*faceV[s])[0])[0]],
lagrange[(*(*faceV[s])[0])[0]] +
lagrange[(*(*faceV[s])[0])[1]]);
u[s][l] = intLegendre[l].compose(lagrange[faceIdx[s][0][1]] -
lagrange[faceIdx[s][0][0]],
lagrange[faceIdx[s][0][0]] +
lagrange[faceIdx[s][0][1]]);
v[s][l] = legendre[l].compose(p[s]);
v[s][l].mul(lagrange[(*(*faceV[s])[0])[2]]);
v[s][l].mul(lagrange[faceIdx[s][0][2]]);
}
// Differences of grad(u) and grad(v)
vector<Polynomial> gradL1 = lagrange[(*(*faceV[s])[0])[0]].gradient();
vector<Polynomial> gradL2 = lagrange[(*(*faceV[s])[0])[1]].gradient();
vector<Polynomial> gradL1 = lagrange[faceIdx[s][0][0]].gradient();
vector<Polynomial> gradL2 = lagrange[faceIdx[s][0][1]].gradient();
vector<Polynomial> l2GradL1(gradL1);
l2GradL1[0].mul(lagrange[(*(*faceV[s])[0])[1]]);
l2GradL1[1].mul(lagrange[(*(*faceV[s])[0])[1]]);
l2GradL1[2].mul(lagrange[(*(*faceV[s])[0])[1]]);
l2GradL1[0].mul(lagrange[faceIdx[s][0][1]]);
l2GradL1[1].mul(lagrange[faceIdx[s][0][1]]);
l2GradL1[2].mul(lagrange[faceIdx[s][0][1]]);
vector<Polynomial> l1GradL2(gradL2);
l1GradL2[0].mul(lagrange[(*(*faceV[s])[0])[0]]);
l1GradL2[1].mul(lagrange[(*(*faceV[s])[0])[0]]);
l1GradL2[2].mul(lagrange[(*(*faceV[s])[0])[0]]);
l1GradL2[0].mul(lagrange[faceIdx[s][0][0]]);
l1GradL2[1].mul(lagrange[faceIdx[s][0][0]]);
l1GradL2[2].mul(lagrange[faceIdx[s][0][0]]);
subGrad[s] = new vector<Polynomial>(3);
(*subGrad[s])[0].sub(l1GradL2[0]);
......@@ -155,11 +153,11 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
}
// Face Basis
for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nEdge;
for(size_t s = 0; s < nRefSpace; s++){
size_t i = nEdge;
// Type 1
for(unsigned int l1 = 1; l1 < order; l1++){
for(size_t l1 = 1; l1 < order; l1++){
for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
vector<Polynomial> tmp1 = v[s][l2].gradient();
vector<Polynomial> tmp2 = u[s][l1].gradient();
......@@ -183,7 +181,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
}
// Type 2
for(unsigned int l1 = 1; l1 < order; l1++){
for(size_t l1 = 1; l1 < order; l1++){
for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
vector<Polynomial> tmp1 = v[s][l2].gradient();
vector<Polynomial> tmp2 = u[s][l1].gradient();
......@@ -221,13 +219,13 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
}
// Free Temporary Sapce //
for(unsigned int s = 0; s < nRefSpace; s++)
for(size_t s = 0; s < nRefSpace; s++)
delete[] u[s];
for(unsigned int s = 0; s < nRefSpace; s++)
for(size_t s = 0; s < nRefSpace; s++)
delete[] v[s];
for(unsigned int s = 0; s < nRefSpace; s++)
for(size_t s = 0; s < nRefSpace; s++)
delete subGrad[s];
delete[] legendre;
......@@ -236,22 +234,19 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
delete[] u;
delete[] v;
delete[] subGrad;
*/
}
TriEdgeBasis::~TriEdgeBasis(void){
/*
// ReferenceSpace //
delete refSpace;
// Basis //
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
for(size_t i = 0; i < nRefSpace; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
delete[] basis[i];
}
delete[] basis;
*/
}
FunctionSpace/TriEdgeBasis.h
+ 1
− 1
View file @ f529232c
......@@ -20,7 +20,7 @@ class TriEdgeBasis: public BasisHierarchical1From{
//! @param order The order of the Basis
//!
//! Returns a new Edge-Basis for Triangles of the given order
TriEdgeBasis(unsigned int order);
TriEdgeBasis(size_t order);
//! Deletes this Basis
//!
......
FunctionSpace/TriNedelecBasis.cpp
+ 15
− 20
View file @ f529232c
......@@ -4,13 +4,12 @@
using namespace std;
TriNedelecBasis::TriNedelecBasis(void){
/*
// Reference Space //
refSpace = new TriReferenceSpace;
nRefSpace = refSpace->getNReferenceSpace();
const vector<const vector<const vector<unsigned int>*>*>&
edgeV = refSpace->getAllEdge();
const vector<vector<vector<size_t> > >&
edgeIdx = refSpace->getEdgeNodeIndex();
// Set Basis Type //
order = 0;
......@@ -39,23 +38,22 @@ TriNedelecBasis::TriNedelecBasis(void){
// Basis //
basis = new vector<Polynomial>**[nRefSpace];
for(unsigned int s = 0; s < nRefSpace; s++)
for(size_t s = 0; s < nRefSpace; s++)
basis[s] = new vector<Polynomial>*[nFunction];
// Edge Based (Nedelec) //
for(unsigned int s = 0; s < nRefSpace; s++){
for(unsigned int e = 0; e < 3; e++){
vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].gradient();
for(size_t s = 0; s < nRefSpace; s++){
for(size_t e = 0; e < 3; e++){
vector<Polynomial> tmp1 = lagrange[edgeIdx[s][e][1]].gradient();
vector<Polynomial> tmp2 = lagrange[edgeIdx[s][e][0]].gradient();
tmp1[0].mul(lagrange[(*(*edgeV[s])[e])[0]]);
tmp1[1].mul(lagrange[(*(*edgeV[s])[e])[0]]);
tmp1[2].mul(lagrange[(*(*edgeV[s])[e])[0]]);
tmp1[0].mul(lagrange[edgeIdx[s][e][0]]);
tmp1[1].mul(lagrange[edgeIdx[s][e][0]]);
tmp1[2].mul(lagrange[edgeIdx[s][e][0]]);
tmp2[0].mul(lagrange[(*(*edgeV[s])[e])[1]]);
tmp2[1].mul(lagrange[(*(*edgeV[s])[e])[1]]);
tmp2[2].mul(lagrange[(*(*edgeV[s])[e])[1]]);
tmp2[0].mul(lagrange[edgeIdx[s][e][1]]);
tmp2[1].mul(lagrange[edgeIdx[s][e][1]]);
tmp2[2].mul(lagrange[edgeIdx[s][e][1]]);
tmp2[0].sub(tmp1[0]);
tmp2[1].sub(tmp1[1]);
......@@ -64,22 +62,19 @@ TriNedelecBasis::TriNedelecBasis(void){
basis[s][e] = new vector<Polynomial>(tmp2);
}
}
*/
}
TriNedelecBasis::~TriNedelecBasis(void){
/*
// ReferenceSpace //
delete refSpace;
// Basis //
for(unsigned int i = 0; i < nRefSpace; i++){
for(unsigned int j = 0; j < nFunction; j++)
for(size_t i = 0; i < nRefSpace; i++){
for(size_t j = 0; j < nFunction; j++)
delete basis[i][j];
delete[] basis[i];
}
delete[] basis;
*/
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment