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

Bug Fix: ReferenceSpace compatible with Gmsh MElements (automatic)

parent 8be32242
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
FunctionSpace/TetNodeBasis.cpp
+ 12
− 16
View file @ ac0a46df
...@@ -28,9 +28,9 @@ TetNodeBasis::TetNodeBasis(unsigned int order){ ...@@ -28,9 +28,9 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
nFunction = nVertex + nEdge + nFace + nCell; nFunction = nVertex + nEdge + nFace + nCell;
// Alloc Temporary Space // // Alloc Temporary Space //
const unsigned int orderMinus = order - 1; const int orderMinus = order - 1;
const unsigned int orderMinusTwo = order - 2; const int orderMinusTwo = order - 2;
const unsigned int orderMinusThree = order - 3; const int orderMinusThree = order - 3;
Polynomial* legendre = new Polynomial[order]; Polynomial* legendre = new Polynomial[order];
Polynomial* sclLegendre = new Polynomial[order]; Polynomial* sclLegendre = new Polynomial[order];
...@@ -72,11 +72,11 @@ TetNodeBasis::TetNodeBasis(unsigned int order){ ...@@ -72,11 +72,11 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
} }
// Edge Based // // Edge Based //
for(unsigned int s = 0; s < 2; s++){ for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nVertex; unsigned int i = nVertex;
for(unsigned int l = 1; l < order; l++){ for(unsigned int l = 1; l < order; l++){
for(unsigned int e = 0; e < 6; e++){ for(int e = 0; e < 6; e++){
(*(*basis)[s])[i] = (*(*basis)[s])[i] =
new Polynomial(intLegendre[l].compose new Polynomial(intLegendre[l].compose
(lagrange[(*(*edgeV[s])[e])[0]] - (lagrange[(*(*edgeV[s])[e])[0]] -
...@@ -84,19 +84,18 @@ TetNodeBasis::TetNodeBasis(unsigned int order){ ...@@ -84,19 +84,18 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
, ,
lagrange[(*(*edgeV[s])[e])[0]] + lagrange[(*(*edgeV[s])[e])[0]] +
lagrange[(*(*edgeV[s])[e])[1]])); lagrange[(*(*edgeV[s])[e])[1]]));
i++; i++;
} }
} }
} }
// Face Based // // Face Based //
for(unsigned int s = 0; s < 6; s++){ for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nVertex + nEdge; unsigned int i = nVertex + nEdge;
for(unsigned int l1 = 1; l1 < orderMinus; l1++){ for(int l1 = 1; l1 < orderMinus; l1++){
for(unsigned int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){ for(int l2 = 0; l1 + l2 - 1 < orderMinusTwo; l2++){
for(unsigned int f = 0; f < 4; f++){ for(int f = 0; f < 4; f++){
Polynomial sum = Polynomial sum =
lagrange[(*(*faceV[s])[f])[0]] + lagrange[(*(*faceV[s])[f])[0]] +
lagrange[(*(*faceV[s])[f])[1]] + lagrange[(*(*faceV[s])[f])[1]] +
...@@ -116,7 +115,6 @@ TetNodeBasis::TetNodeBasis(unsigned int order){ ...@@ -116,7 +115,6 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
* *
sclLegendre[l2].compose sclLegendre[l2].compose
(lagrange[(*(*faceV[s])[f])[2]] * 2 - sum, sum)); (lagrange[(*(*faceV[s])[f])[2]] * 2 - sum, sum));
i++; i++;
} }
} }
...@@ -134,10 +132,9 @@ TetNodeBasis::TetNodeBasis(unsigned int order){ ...@@ -134,10 +132,9 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
for(unsigned int s = 0; s < nRefSpace; s++){ for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nVertex + nEdge + nFace; unsigned int i = nVertex + nEdge + nFace;
for(unsigned int l1 = 1; l1 < orderMinusTwo; l1++){ for(int l1 = 1; l1 < orderMinusTwo; l1++){
for(unsigned int l2 = 0; l2 + l1 - 1 < orderMinusThree; l2++){ for(int l2 = 0; l2 + l1 - 1 < orderMinusThree; l2++){
for(unsigned int l3 = 0; l3 + l2 + l1 - 1 < orderMinusThree; l3++){ for(int l3 = 0; l3 + l2 + l1 - 1 < orderMinusThree; l3++){
(*(*basis)[s])[i] = (*(*basis)[s])[i] =
new Polynomial(intLegendre[l1].compose(sub, add) * new Polynomial(intLegendre[l1].compose(sub, add) *
lagrange[2] * lagrange[2] *
...@@ -145,7 +142,6 @@ TetNodeBasis::TetNodeBasis(unsigned int order){ ...@@ -145,7 +142,6 @@ TetNodeBasis::TetNodeBasis(unsigned int order){
oneMinusFour) * oneMinusFour) *
lagrange[3] * lagrange[3] *
legendre[l3].compose(twoFourMinusOne)); legendre[l3].compose(twoFourMinusOne));
i++; i++;
} }
} }
......
FunctionSpace/TetReferenceSpace.cpp
+ 13
− 17
View file @ ac0a46df
#include <sstream> #include <sstream>
#include "TetReferenceSpace.h" #include "TetReferenceSpace.h"
#include "MTetrahedron.h"
using namespace std; using namespace std;
...@@ -11,28 +12,23 @@ TetReferenceSpace::TetReferenceSpace(void){ ...@@ -11,28 +12,23 @@ TetReferenceSpace::TetReferenceSpace(void){
nEdge = 6; nEdge = 6;
refEdge = new unsigned int*[nEdge]; refEdge = new unsigned int*[nEdge];
for(unsigned int i = 0; i < nEdge; i++) for(unsigned int i = 0; i < nEdge; i++){
refEdge[i] = new unsigned int[2]; refEdge[i] = new unsigned int[2];
refEdge[i][0] = MTetrahedron::edges_tetra(i, 0);
refEdge[0][0] = 0; refEdge[0][1] = 1; refEdge[i][1] = MTetrahedron::edges_tetra(i, 1);
refEdge[1][0] = 1; refEdge[1][1] = 2; }
refEdge[2][0] = 2; refEdge[2][1] = 0;
refEdge[3][0] = 0; refEdge[3][1] = 3;
refEdge[4][0] = 1; refEdge[4][1] = 3;
refEdge[5][0] = 2; refEdge[5][1] = 3;
// Face Definition // // Face Definition //
nFace = 4; nFace = 4;
refFace = new unsigned int*[nFace]; refFace = new unsigned int*[nFace];
for(unsigned int i = 0; i < nFace; i++) for(unsigned int i = 0; i < nFace; i++){
refFace[i] = new unsigned int[3]; refFace[i] = new unsigned int[3];
refFace[i][0] = MTetrahedron::faces_tetra(i, 0);
refFace[i][1] = MTetrahedron::faces_tetra(i, 1);
refFace[i][2] = MTetrahedron::faces_tetra(i, 2);
}
refFace[0][0] = 0; refFace[0][1] = 1; refFace[0][2] = 2;
refFace[1][0] = 0; refFace[1][1] = 1; refFace[1][2] = 3;
refFace[2][0] = 1; refFace[2][1] = 2; refFace[2][2] = 3;
refFace[3][0] = 2; refFace[3][1] = 0; refFace[3][2] = 3;
// Init All // // Init All //
init(); init();
} }
......
FunctionSpace/TriEdgeBasis.cpp
+ 9
− 9
View file @ ac0a46df
...@@ -25,8 +25,8 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){ ...@@ -25,8 +25,8 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
nFunction = nVertex + nEdge + nFace + nCell; nFunction = nVertex + nEdge + nFace + nCell;
// Alloc Some Space // // Alloc Some Space //
const unsigned int orderPlus = order + 1; const int orderPlus = order + 1;
const unsigned int orderMinus = order - 1; const int orderMinus = order - 1;
Polynomial* legendre = new Polynomial[orderPlus]; Polynomial* legendre = new Polynomial[orderPlus];
Polynomial* intLegendre = new Polynomial[orderPlus]; Polynomial* intLegendre = new Polynomial[orderPlus];
...@@ -57,7 +57,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){ ...@@ -57,7 +57,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
// Edge Based (Nedelec) // // Edge Based (Nedelec) //
for(unsigned int s = 0; s < nRefSpace; s++){ for(unsigned int s = 0; s < nRefSpace; s++){
for(unsigned int e = 0; e < 3; e++){ for(int e = 0; e < 3; e++){
vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient(); vector<Polynomial> tmp1 = lagrange[(*(*edgeV[s])[e])[1]].gradient();
vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].gradient(); vector<Polynomial> tmp2 = lagrange[(*(*edgeV[s])[e])[0]].gradient();
...@@ -82,8 +82,8 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){ ...@@ -82,8 +82,8 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
for(unsigned int s = 0; s < nRefSpace; s++){ for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = 3; unsigned int i = 3;
for(unsigned int l = 1; l < orderPlus; l++){ for(int l = 1; l < orderPlus; l++){
for(unsigned int e = 0; e < 3; e++){ for(int e = 0; e < 3; e++){
(*(*basis)[s])[i] = (*(*basis)[s])[i] =
new vector<Polynomial> new vector<Polynomial>
((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] - ((intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[0]] -
...@@ -101,7 +101,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){ ...@@ -101,7 +101,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
// Preliminaries // Preliminaries
const Polynomial p = lagrange[2] * 2 - Polynomial(1, 0, 0, 0); const Polynomial p = lagrange[2] * 2 - Polynomial(1, 0, 0, 0);
for(unsigned int l = 0; l < orderPlus; l++){ for(int l = 0; l < orderPlus; l++){
u[l] = intLegendre[l].compose(lagrange[1] - lagrange[0], u[l] = intLegendre[l].compose(lagrange[1] - lagrange[0],
lagrange[0] + lagrange[1]); lagrange[0] + lagrange[1]);
v[l] = legendre[l].compose(p); v[l] = legendre[l].compose(p);
...@@ -131,7 +131,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){ ...@@ -131,7 +131,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
// Type 1 // Type 1
for(unsigned int l1 = 1; l1 < order; l1++){ for(unsigned int l1 = 1; l1 < order; l1++){
for(unsigned int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){ for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
vector<Polynomial> tmp1 = v[l2].gradient(); vector<Polynomial> tmp1 = v[l2].gradient();
vector<Polynomial> tmp2 = u[l1].gradient(); vector<Polynomial> tmp2 = u[l1].gradient();
...@@ -155,7 +155,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){ ...@@ -155,7 +155,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
// Type 2 // Type 2
for(unsigned int l1 = 1; l1 < order; l1++){ for(unsigned int l1 = 1; l1 < order; l1++){
for(unsigned int l2 = 0; l2 + l1 - 1 < orderMinus; l2++){ for(int l2 = 0; l2 + (int)l1 - 1 < orderMinus; l2++){
vector<Polynomial> tmp1 = v[l2].gradient(); vector<Polynomial> tmp1 = v[l2].gradient();
vector<Polynomial> tmp2 = u[l1].gradient(); vector<Polynomial> tmp2 = u[l1].gradient();
...@@ -178,7 +178,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){ ...@@ -178,7 +178,7 @@ TriEdgeBasis::TriEdgeBasis(unsigned int order){
} }
// Type 3 // Type 3
for(unsigned int l = 0; l < orderMinus; l++){ for(int l = 0; l < orderMinus; l++){
vector<Polynomial> subGradL1L2V(subGradL1L2); vector<Polynomial> subGradL1L2V(subGradL1L2);
subGradL1L2V[0].mul(v[l]); subGradL1L2V[0].mul(v[l]);
......
FunctionSpace/TriNodeBasis.cpp
+ 7
− 6
View file @ ac0a46df
...@@ -25,7 +25,7 @@ TriNodeBasis::TriNodeBasis(unsigned int order){ ...@@ -25,7 +25,7 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
nFunction = nVertex + nEdge + nFace + nCell; nFunction = nVertex + nEdge + nFace + nCell;
// Alloc Some Space // // Alloc Some Space //
const unsigned int orderMinus = order - 1; const int orderMinus = order - 1;
Polynomial* legendre = new Polynomial[order]; Polynomial* legendre = new Polynomial[order];
Polynomial* intLegendre = new Polynomial[order]; Polynomial* intLegendre = new Polynomial[order];
...@@ -64,7 +64,7 @@ TriNodeBasis::TriNodeBasis(unsigned int order){ ...@@ -64,7 +64,7 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
unsigned int i = nVertex; unsigned int i = nVertex;
for(unsigned int l = 1; l < order; l++){ for(unsigned int l = 1; l < order; l++){
for(unsigned int e = 0; e < 3; e++){ for(int e = 0; e < 3; e++){
(*(*basis)[s])[i] = (*(*basis)[s])[i] =
new Polynomial(intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[1]] - new Polynomial(intLegendre[l].compose(lagrange[(*(*edgeV[s])[e])[1]] -
lagrange[(*(*edgeV[s])[e])[0]] lagrange[(*(*edgeV[s])[e])[0]]
...@@ -78,16 +78,17 @@ TriNodeBasis::TriNodeBasis(unsigned int order){ ...@@ -78,16 +78,17 @@ TriNodeBasis::TriNodeBasis(unsigned int order){
// Cell Based // // Cell Based //
const Polynomial p = (lagrange[2] * 2) - Polynomial(1, 0, 0, 0); const Polynomial p = (lagrange[2] * 2) - Polynomial(1, 0, 0, 0);
const unsigned int orderMinusTwo = order - 2; const int orderMinusTwo = order - 2;
for(unsigned int s = 0; s < nRefSpace; s++){ for(unsigned int s = 0; s < nRefSpace; s++){
unsigned int i = nVertex + nEdge; unsigned int i = nVertex + nEdge;
for(unsigned int l1 = 1; l1 < orderMinus; l1++){ for(int l1 = 1; l1 < orderMinus; l1++){
for(unsigned int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){ for(int l2 = 0; l2 + l1 - 1 < orderMinusTwo; l2++){
(*(*basis)[s])[i] = (*(*basis)[s])[i] =
new Polynomial(intLegendre[l1].compose(lagrange[1] - lagrange[0], new Polynomial(intLegendre[l1].compose(lagrange[1] - lagrange[0],
lagrange[1] + lagrange[0]) * lagrange[1] + lagrange[0])
*
legendre[l2].compose(p) * lagrange[2]); legendre[l2].compose(p) * lagrange[2]);
i++; i++;
......
FunctionSpace/TriReferenceSpace.cpp
+ 6
− 6
View file @ ac0a46df
#include <sstream> #include <sstream>
#include "TriReferenceSpace.h" #include "TriReferenceSpace.h"
#include "MTriangle.h"
using namespace std; using namespace std;
...@@ -11,12 +12,11 @@ TriReferenceSpace::TriReferenceSpace(void){ ...@@ -11,12 +12,11 @@ TriReferenceSpace::TriReferenceSpace(void){
nEdge = 3; nEdge = 3;
refEdge = new unsigned int*[nEdge]; refEdge = new unsigned int*[nEdge];
for(unsigned int i = 0; i < nEdge; i++) for(unsigned int i = 0; i < nEdge; i++){
refEdge[i] = new unsigned int[2]; refEdge[i] = new unsigned int[2];
refEdge[i][0] = MTriangle::edges_tri(i, 0);
refEdge[0][0] = 0; refEdge[0][1] = 1; refEdge[i][1] = MTriangle::edges_tri(i, 1);
refEdge[1][0] = 1; refEdge[1][1] = 2; }
refEdge[2][0] = 2; refEdge[2][1] = 0;
// Face Definition // // Face Definition //
nFace = 0; nFace = 0;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment