From 01ca5b4caafa447c4cc65b6092737895621c2ce3 Mon Sep 17 00:00:00 2001
From: Tuomas Karna <tuomas.karna@uclouvain.be>
Date: Fri, 5 Mar 2010 17:02:52 +0000
Subject: [PATCH] dg: added 2nd order Prism18
---
Geo/GModel.cpp | 2 +-
Geo/MPrism.cpp | 4 +-
Geo/MPrism.h | 9 +--
Numeric/GaussQuadraturePri.cpp | 2 +-
Numeric/polynomialBasis.cpp | 105 ++++++++++++++++++++-----------
Solver/TESTCASES/Diffusion3D.lua | 23 +++----
Solver/TESTCASES/mixed.geo | 6 +-
Solver/dgDofContainer.cpp | 4 ++
8 files changed, 92 insertions(+), 63 deletions(-)
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 5aff9d50c4..346edff76f 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1334,7 +1334,7 @@ void GModel::load(std::string fileName)
{
GModel *temp = GModel::current();
GModel::setCurrent(this);
- MergeFile(fileName.c_str());
+ MergeFile(fileName.c_str(),true);
GModel::setCurrent(temp);
}
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index dbe7e9b8ac..4728e9d722 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -37,14 +37,14 @@ const polynomialBasis* MPrism::getFunctionSpace(int o) const
switch (order) {
case 1: return &polynomialBases::find(MSH_PRI_6);
case 2: return &polynomialBases::find(MSH_PRI_18);
- default: Msg::Error("Order %d tetrahedron function space not implemented", order);
+ default: Msg::Error("Order %d prism function space not implemented", order);
}
}
else {
switch (order) {
case 1: return &polynomialBases::find(MSH_PRI_6);
case 2: return &polynomialBases::find(MSH_PRI_18);
- default: Msg::Error("Order %d tetrahedron function space not implemented", order);
+ default: Msg::Error("Order %d prism function space not implemented", order);
}
}
return 0;
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index 3d1220f472..97a796eb18 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -131,7 +131,7 @@ class MPrism : public MElement {
}
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
virtual int getVolumeSign();
- virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
+/* virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
{
s[0] = (1. - u - v) * (1. - w) * 0.5;
s[1] = u * (1. - w) * 0.5;
@@ -139,9 +139,10 @@ class MPrism : public MElement {
s[3] = (1. - u - v) * (1. + w) * 0.5;
s[4] = u * (1. + w) * 0.5;
s[5] = v * (1. + w) * 0.5;
- }
- virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+ }*/
+/* virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
{
+ printf("here\n");
s[0][0] = -0.5 * (1. - w) ;
s[0][1] = -0.5 * (1. - w) ;
s[0][2] = -0.5 * (1. - u - v);
@@ -160,7 +161,7 @@ class MPrism : public MElement {
s[5][0] = 0. ;
s[5][1] = 0.5 * (1. + w) ;
s[5][2] = 0.5 * v ;
- }
+ }*/
virtual bool isInside(double u, double v, double w)
{
double tol = _isInsideTolerance;
diff --git a/Numeric/GaussQuadraturePri.cpp b/Numeric/GaussQuadraturePri.cpp
index 3996d29543..7696549480 100644
--- a/Numeric/GaussQuadraturePri.cpp
+++ b/Numeric/GaussQuadraturePri.cpp
@@ -10,7 +10,7 @@
IntPt *getGQPriPts(int order);
int getNGQPriPts(int order);
-IntPt * GQP[20] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
+IntPt * GQP[28] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
IntPt *getGQPriPts(int order)
{
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index d63c948e43..b2a59a6690 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -244,17 +244,15 @@ static fullMatrix<double> generatePascalPrism(int order)
int index = 0;
fullMatrix<double> lineMonoms = generate1DMonomials(order);
fullMatrix<double> triMonoms = generatePascalTriangle(order);
-// printf("nb: %d nT: %d %d nL: %d %d\n",nbMonomials,triMonoms.size1(),(order+1)*(order+2)/2,lineMonoms.size1(),order+1);
for (int j = 0; j < lineMonoms.size1(); j++) {
for (int i = 0; i < triMonoms.size1(); i++) {
- monomials(index,0) = triMonoms(i,0);
- monomials(index,1) = triMonoms(i,1);
- monomials(index,2) = lineMonoms(j,0);
- index ++;
-// printf("%d: %f %f %f\n",index,triMonoms(i,0),triMonoms(i,1),lineMonoms(j,0));
+ monomials(index,0) = triMonoms(i,0);
+ monomials(index,1) = triMonoms(i,1);
+ monomials(index,2) = lineMonoms(j,0);
+ index ++;
}
}
-// monomials.print();
+ monomials.print("Pri monoms");
return monomials;
}
@@ -625,26 +623,50 @@ static fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
static fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip)
{
+ const double prism18Pts[18][3] = {
+ {0, 0, -1}, // 0
+ {1, 0, -1}, // 1
+ {0, 1, -1}, // 2
+ {0, 0, 1}, // 3
+ {1, 0, 1}, // 4
+ {0, 1, 1}, // 5
+ {0.5, 0, -1}, // 6
+ {0, 0.5, -1}, // 7
+ {0, 0, 0}, // 8
+ {0.5, 0.5, -1}, // 9
+ {1, 0, 0}, // 10
+ {0, 1, 0}, // 11
+ {0.5, 0, 1}, // 12
+ {0, 0.5, 1}, // 13
+ {0.5, 0.5, 1}, // 14
+ {0.5, 0, 0}, // 15
+ {0, 0.5, 0}, // 16
+ {0.5, 0.5, 0}, // 17
+ };
+
int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
fullMatrix<double> point(nbPoints, 3);
- double overOrder = (order == 0 ? 1. : 1. / order);
int index = 0;
fullMatrix<double> triPoints = gmshGeneratePointsTriangle(order,false);
fullMatrix<double> linePoints = generate1DPoints(order);
-// printf("nb: %d nT: %d %d nL: %d %d\n",nbPoints,triPoints.size1(),(order+1)*(order+2)/2,linePoints.size1(),order+1);
- for (int j = 0; j < linePoints.size1(); j++) {
- for (int i = 0; i < triPoints.size1(); i++) {
- point(index,0) = triPoints(i,0);
- point(index,1) = triPoints(i,1);
- point(index,2) = linePoints(j,0);
- index ++;
-// printf("%d: %f %f %f\n",index,triPoints(i,0),triPoints(i,1),linePoints(j,0));
+
+ if (order == 2)
+ for (int i =0; i<18; i++)
+ for (int j=0; j<3;j++)
+ point(i,j) = prism18Pts[i][j];
+ else
+ for (int j = 0; j <linePoints.size1() ; j++) {
+ for (int i = 0; i < triPoints.size1(); i++) {
+ point(index,0) = triPoints(i,0);
+ point(index,1) = triPoints(i,1);
+ point(index,2) = linePoints(j,0);
+ index ++;
+ }
}
- }
-// point.print();
+
+// point.print("Pri ipts");
- point.scale(overOrder);
return point;
}
@@ -798,26 +820,31 @@ static void generate3dFaceClosure(polynomialBasis::clCont &closure, int order)
static void getFaceClosurePrism(int iFace, int iSign, int iRotate, std::vector<int> &closure,
int order)
{
- if (order > 1)
+ if (order > 2)
Msg::Error("FaceClosure not implemented for prisms of order %d",order);
bool isTriangle = iFace<2;
int nNodes = isTriangle ? (order+1)*(order+2)/2 : (order+1)*(order+1);
closure.clear();
closure.resize(nNodes);
- switch (order){
- case 0:
+ if (order==0) {
closure[0] = 0;
- break;
- default:
-// int face[4][3] = {{-3, -2, -1}, {1, -6, 4}, {-4, 5, 3}, {6, 2, -5}};
- int order1node[5][4] = {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {0, 3, 5, 2}, {1, 2, 5, 4}};
- int nVertex = isTriangle ? 3 : 4;
+ return;
+ }
+ int order1node[5][4] = {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {0, 3, 5, 2}, {1, 2, 5, 4}};
+ int order2node[5][5] = {{7, 9, 6, -1, -1}, {12, 14, 13, -1, -1}, {6, 10, 12, 8, 15}, {8, 13, 11, 7, 16}, {9, 11, 14, 10, 17}};
+// int order2node[5][4] = {{7, 9, 6, -1}, {12, 14, 13, -1}, {6, 10, 12, 8}, {8, 13, 11, 7}, {9, 11, 14, 10}};
+ int nVertex = isTriangle ? 3 : 4;
+ for (int i = 0; i < nVertex; ++i){
+ int k = (nVertex + (iSign * i) + iRotate) % nVertex; //- iSign * iRotate
+ closure[i] = order1node[iFace][k];
+ }
+ if (order==2) {
for (int i = 0; i < nVertex; ++i){
- int k;
- k = (nVertex + (iSign * i) + iRotate) % nVertex; //- iSign * iRotate
- closure[i] = order1node[iFace][k];
+ int k = (nVertex + (iSign==-1?-1:0) + (iSign * i) + iRotate) % nVertex; //- iSign * iRotate
+ closure[nVertex+i] = order2node[iFace][k];
}
- break;
+ if (!isTriangle)
+ closure[nNodes-1] = order2node[iFace][4]; // center
}
}
@@ -828,11 +855,9 @@ static void generate3dFaceClosurePrism(polynomialBasis::clCont &closure, int ord
for (int iRotate = 0; iRotate < 4; iRotate++){
for (int iSign = 1; iSign >= -1; iSign -= 2){
for (int iFace = 0; iFace < 5; iFace++){
- std::vector<int> closure_face;
- getFaceClosurePrism(iFace, iSign, iRotate, closure_face, order);
- closure.push_back(closure_face);
-// for(int i=0;i < closure_face.size(); i++) { printf("%d ",closure_face.at(i)); }
-// printf("\n");
+ std::vector<int> closure_face;
+ getFaceClosurePrism(iFace, iSign, iRotate, closure_face, order);
+ closure.push_back(closure_face);
}
}
}
@@ -1075,6 +1100,14 @@ const polynomialBasis &polynomialBases::find(int tag)
break;
}
F.coefficients = generateLagrangeMonomialCoefficients(F.monomials, F.points);
+// printf("Case: %d coeffs:\n",tag);
+// for (int i = 0; i<F.coefficients.size1(); i++) {
+// for (int j = 0; j<F.coefficients.size2(); j++) {
+// printf("%4.1f ",F.coefficients(i,j));
+// }
+// printf("\n");
+// }
+
fs.insert(std::make_pair(tag, F));
return fs[tag];
}
diff --git a/Solver/TESTCASES/Diffusion3D.lua b/Solver/TESTCASES/Diffusion3D.lua
index d53fb73abb..ac7550ef65 100644
--- a/Solver/TESTCASES/Diffusion3D.lua
+++ b/Solver/TESTCASES/Diffusion3D.lua
@@ -6,7 +6,7 @@ model = GModel()
model:load('mixed.geo')
model:load('mixed.msh')
dg = dgSystemOfEquations(model)
-dg:setOrder(1)
+dg:setOrder(2)
-- initial condition
function initial_condition( xyz , f )
@@ -37,17 +37,6 @@ dg:setConservationLaw(law)
-- boundary condition
outside=fullMatrix(1,1)
outside:set(0,0,0.15)
--- freestream=law:newOutsideValueBoundary(functionLua(1,'initial_condition',{'XYZ'}):getName())
-
--- law:addBoundaryCondition('wall',law:new0FluxBoundary())
--- law:addBoundaryCondition('inlet',freestream)
--- law:addBoundaryCondition('outlet',law:newSymmetryBoundary())
--- law:addBoundaryCondition('symmetry',law:newSymmetryBoundary())
-
--- law:addBoundaryCondition('wall',law:new0FluxBoundary())
--- law:addBoundaryCondition('inlet',law:new0FluxBoundary())
--- law:addBoundaryCondition('outlet',law:new0FluxBoundary())
--- law:addBoundaryCondition('symmetry',law:new0FluxBoundary())
law:addBoundaryCondition('side1',law:new0FluxBoundary())
law:addBoundaryCondition('side2',law:new0FluxBoundary())
@@ -60,10 +49,10 @@ dg:setup()
dg:L2Projection(functionLua(1,'initial_condition',{'XYZ'}):getName())
-dt = 0.4
+dt = 0.0001
--CFL = 20 -- good for lc=0.1 mesh
-CFL = 10
-dt = CFL*dg:computeInvSpectralRadius()
+CFL = 1
+-- dt = CFL*dg:computeInvSpectralRadius()
-- export_interval = 5*dt
export_interval = 0.05
end_time = 3
@@ -78,11 +67,13 @@ export_count = 0
io.write(string.format('Export: %d iter: %4d t: %.6f n: %.8f dt:%.8f\n',export_count,0,0,0,dt))
dg:exportSolution(output_pattern .. string.format("%05d",export_count))
for i=1,10000 do
- dt = CFL*dg:computeInvSpectralRadius()
+-- dt = CFL*dg:computeInvSpectralRadius()
norm = dg:RK44(dt)
time = i*dt
io.write('.')
io.flush()
+-- io.write('\n')
+-- io.write(string.format('Export: %d iter: %4d t: %.6f n: %.8f dt:%.8f\n',export_count,i,time,norm,dt))
if (time >= next_export_time) then
export_count = export_count +1
io.write('\n')
diff --git a/Solver/TESTCASES/mixed.geo b/Solver/TESTCASES/mixed.geo
index f5c0427daf..c129eea5fb 100644
--- a/Solver/TESTCASES/mixed.geo
+++ b/Solver/TESTCASES/mixed.geo
@@ -2,7 +2,7 @@
C = 1;
Lup = 1;
L = 1.;
-lc = .3;
+lc = 0.3;
Point(1) = {0.0, 0.0, -Lup, lc};
Point(2) = {C , 0.0, -Lup, lc};
@@ -23,7 +23,7 @@ outtet[] = Extrude {0,0,0.7*L} { Surface{6};};
// Printf("volume = %g", outtet[1]);
// Printf("side surfaces = %g %g %g %g", outtet[2], outtet[3], outtet[4], outtet[5]);
-outpri[]= Extrude {0,0,0.5*L}{ Surface{outtet[0]}; Layers{4};Recombine;};
+outpri[]= Extrude {0,0,0.5*L}{ Surface{outtet[0]}; Layers{Ceil(0.5*L/lc)};Recombine;};
// outv[]= Extrude {0,0,0.5*L}{ Surface{7}; Layers{5};Recombine;};
// Printf("top surface = %g", outpri[0]);
// Printf("volume = %g", outpri[1]);
@@ -32,7 +32,7 @@ outpri[]= Extrude {0,0,0.5*L}{ Surface{outtet[0]}; Layers{4};Recombine;};
Mesh.Algorithm3D=4; // frontal [lobes]
Physical Surface("top") = {outpri[0]};
-Physical Surface("bottom") = {5};
+Physical Surface("bottom") = {6};
Physical Surface("side1") = {outpri[2],outtet[2]};
Physical Surface("side2") = {outpri[3],outtet[3]};
Physical Surface("side3") = {outpri[4],outtet[4]};
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
index 271b347ad6..b976dd25f3 100644
--- a/Solver/dgDofContainer.cpp
+++ b/Solver/dgDofContainer.cpp
@@ -304,6 +304,10 @@ void dgDofContainer::exportMsh(const std::string name)
if(Msg::GetCommSize()>1)
name_oss<<"_"<<Msg::GetCommRank();
FILE *f = fopen (name_oss.str().c_str(),"w");
+ if(!f){
+ Msg::Error("Unable to open export file '%s'", name.c_str());
+ }
+
int COUNT = 0;
for (int i=0;i < _groups.getNbElementGroups() ;i++){
COUNT += _groups.getElementGroup(i)->getNbElements();
--
GitLab