From 81f2ac09fd3080a5c317f295d4e6712fd3f06442 Mon Sep 17 00:00:00 2001
From: Samuel Melchior <samuel.melchior@uclouvain.be>
Date: Tue, 15 Dec 2009 09:13:00 +0000
Subject: [PATCH] 1D lua
---
Solver/TESTCASES/Advection1D.lua | 53 ++++++++++++++++++++++++++++++++
Solver/TESTCASES/WavePulse.lua | 2 +-
Solver/dgGroupOfElements.cpp | 5 ++-
Solver/dgSystemOfEquations.cpp | 4 +--
4 files changed, 60 insertions(+), 4 deletions(-)
create mode 100644 Solver/TESTCASES/Advection1D.lua
diff --git a/Solver/TESTCASES/Advection1D.lua b/Solver/TESTCASES/Advection1D.lua
new file mode 100644
index 0000000000..c85e108ecd
--- /dev/null
+++ b/Solver/TESTCASES/Advection1D.lua
@@ -0,0 +1,53 @@
+model = GModel ()
+model:load ('edge.geo')
+--model:mesh (1)
+--model:save ('edge_new.msh')
+model:load ('edge_new.msh')
+dg = dgSystemOfEquations (model)
+dg:setOrder(1)
+
+
+-- conservation law
+-- advection speed
+v=fullMatrix(3,1);
+v:set(0,0,0.15)
+v:set(1,0,0)
+v:set(2,0,0)
+-- diffusivity
+nu=fullMatrix(1,1);
+nu:set(0,0,0)
+
+dg:setConservationLaw('AdvectionDiffusion',createFunction.constant(v),createFunction.constant(nu))
+
+-- boundary condition
+outside=fullMatrix(1,1)
+outside:set(0,0,0.)
+dg:addBoundaryCondition('Left','OutsideValues',createFunction.constant(outside))
+dg:addBoundaryCondition('Right','OutsideValues',createFunction.constant(outside))
+
+dg:setup()
+
+-- initial condition
+function initial_condition( _x , _f )
+ xyz = fullMatrix(_x)
+ f = fullMatrix(_f)
+ for i=0,xyz:size1()-1 do
+ x = xyz:get(i,0)
+ y = xyz:get(i,1)
+ z = xyz:get(i,2)
+ f:set (i, 0, math.exp(-100*((x+0.8)^2)))
+ end
+end
+dg:L2Projection(createFunction.lua(1,'initial_condition','XYZ'))
+
+dg:exportSolution('output/Adv1D_00000')
+
+-- main loop
+n = 5
+for i=1,100*n do
+ norm = dg:RK44(0.03)
+ if (i % n == 0) then
+ print('iter',i,norm)
+ dg:exportSolution(string.format("output/Adv1D-%05d", i))
+ end
+end
diff --git a/Solver/TESTCASES/WavePulse.lua b/Solver/TESTCASES/WavePulse.lua
index 02dc142988..46da82fa2d 100644
--- a/Solver/TESTCASES/WavePulse.lua
+++ b/Solver/TESTCASES/WavePulse.lua
@@ -49,7 +49,7 @@ for i=1,1000 do
norm = DG:RK44(dt)
print('*** ITER ***',i,norm)
if (i % 10 == 0) then
- DG:exportSolution(string.format("solution-%03d", i))
+ DG:exportSolution(string.format("output/solution-%03d", i))
end
end
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index b2c8a303d6..01059e5a74 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -278,7 +278,7 @@ void dgGroupOfFaces::init(int pOrder) {
}
// there is nothing on the right for boundary groups
if(_fsRight){
- fullMatrix<double> &intRight=_integrationPointsLeft[_closuresIdRight[i]];
+ fullMatrix<double> &intRight=_integrationPointsRight[_closuresIdRight[i]];
for (int j=0; j<intRight.size1(); j++){
_fsRight->df((intRight)(j,0),(intRight)(j,1),(intRight)(j,2),g);
getElementRight(i)->getJacobian ((intRight)(j,0), (intRight)(j,1), (intRight)(j,2), jac);
@@ -314,6 +314,7 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string bo
if(boundaryTag=="")
throw;
_fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
+ _closuresLeft = _fsLeft->vertexClosure;
_fsRight=NULL;
for(int i=0; i<elGroup.getNbElements(); i++){
MElement &el = *elGroup.getElement(i);
@@ -374,6 +375,8 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
switch (_groupLeft.getElement(0)->getDim()) {
case 1 : {
std::map<MVertex*,int> vertexMap;
+ _closuresLeft = _fsLeft->vertexClosure;
+ _closuresRight = _fsRight->vertexClosure;
for(int i=0; i<elGroup.getNbElements(); i++){
MElement &el = *elGroup.getElement(i);
for (int j=0; j<el.getNumVertices(); j++){
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 6372d32dea..8ba9b2b7b1 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -53,7 +53,7 @@ dgSystemOfEquations::dgSystemOfEquations(lua_State *L){
LuaGModel *obj =Luna<LuaGModel>::check(L, 1);
if (!obj)throw;
_gm = obj->getGModel();
- _dimension = _gm->getNumRegions() ? 3 : 2;
+ _dimension = _gm->getNumRegions() ? 3 : _gm->getNumFaces() ? 2 : 1;
_solution = 0;
}
@@ -194,7 +194,7 @@ void dgSystemOfEquations::export_solution_as_is (const std::string &name){
fullMatrix<double> sol = getSolutionProxy (i, iElement);
fprintf(f,"%d %d",num,sol.size1());
for (int k=0;k<sol.size1();++k) {
- fprintf(f,"%12.5E ",sol(k,ICOMP));
+ fprintf(f," %12.5E ",sol(k,ICOMP));
}
fprintf(f,"\n");
}
--
GitLab