Skip to content
Snippets Groups Projects
Commit 81f2ac09 authored by Samuel Melchior's avatar Samuel Melchior
Browse files

1D lua

parent a018dfe3
Branches
Tags
No related merge requests found
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
......@@ -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
......
......@@ -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++){
......
......@@ -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;
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment