dg: Diffusion3D.lua mixed mesh test case

Solver/TESTCASES/Diffusion3D.lua

+model = GModel()
+-- model:load('tetrahedra.geo')
+-- model:load('tetrahedra.msh')
+-- model:load('prisms.geo')
+-- model:load('prisms.msh')
+model:load('mixed.geo')
+model:load('mixed.msh')
+dg = dgSystemOfEquations(model)
+dg:setOrder(1)
+
+-- initial condition
+function initial_condition( xyz , 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(-50*((x-0.4)^2+(y-0.3)^2)))
+--       f:set (i, 0, math.exp(-50*((x-0.4)^2+(z+0.3)^2)))
+   end
+end
+
+-- conservation law
+-- advection speed
+v=fullMatrix(3,1);
+v:set(0,0,0)
+v:set(1,0,0)
+v:set(2,0,0)
+
+nu=fullMatrix(1,1);
+nu:set(0,0,0.01)
+
+-- diffusion function
+
+law = dgConservationLawAdvectionDiffusion(functionConstant(v):getName(),functionConstant(nu):getName())
+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())
+law:addBoundaryCondition('side3',law:new0FluxBoundary())
+law:addBoundaryCondition('side4',law:new0FluxBoundary())
+law:addBoundaryCondition('top',law:newSymmetryBoundary())
+law:addBoundaryCondition('bottom',law:newSymmetryBoundary())
+
+dg:setup()
+
+dg:L2Projection(functionLua(1,'initial_condition',{'XYZ'}):getName())
+
+dt = 0.4
+--CFL = 20 -- good for lc=0.1 mesh
+CFL = 10
+dt = CFL*dg:computeInvSpectralRadius()
+-- export_interval = 5*dt
+export_interval = 0.05
+end_time = 3
+output_dir = 'output'
+output_prefix = 'diff3d'
+output_pattern = output_dir .. '/' .. output_prefix .. '-'
+
+
+-- main loop
+next_export_time = export_interval
+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()
+  norm = dg:RK44(dt)
+  time = i*dt
+  io.write('.')
+  io.flush()
+  if (time >= next_export_time) then
+    export_count = export_count +1
+    io.write('\n')
+    io.write(string.format('Export: %d iter: %4d t: %.6f n: %.8f dt:%.8f\n',export_count,i,time,norm,dt))
+    dg:exportSolution(output_pattern .. string.format("%05d",export_count)) 
+    next_export_time = next_export_time + export_interval
+  end
+  if time >= end_time then
+    break
+  end
+end

Solver/TESTCASES/mixed.geo

+//*********** mixed.geo *************//
+C = 1;
+Lup = 1;
+L = 1.;
+lc = .3;
+
+Point(1) = {0.0, 0.0, -Lup, lc};
+Point(2) = {C  , 0.0, -Lup, lc};
+Point(3) = {C  , C  , -Lup, lc};
+Point(4) = {0.0, C  , -Lup, lc};
+
+Line(1) = {1,2};
+Line(2) = {2,3};
+Line(3) = {3,4};
+Line(4) = {4,1};
+
+
+Line Loop(5) = {1,2,3,4};
+Plane Surface(6) = {5};
+
+outtet[] = Extrude {0,0,0.7*L} { Surface{6};};
+// Printf("top surface = %g", outtet[0]);
+// 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;};
+// outv[]= Extrude {0,0,0.5*L}{ Surface{7}; Layers{5};Recombine;};
+// Printf("top surface = %g", outpri[0]);
+// Printf("volume = %g", outpri[1]);
+// Printf("side surfaces = %g %g %g %g", outpri[2], outpri[3], outpri[4], outpri[5]);
+
+Mesh.Algorithm3D=4; // frontal [lobes]
+
+Physical Surface("top") = {outpri[0]};
+Physical Surface("bottom") = {5};
+Physical Surface("side1") = {outpri[2],outtet[2]};
+Physical Surface("side2") = {outpri[3],outtet[3]};
+Physical Surface("side3") = {outpri[4],outtet[4]};
+Physical Surface("side4") = {outpri[5],outtet[5]};
+Physical Volume("volume") = {outtet[1],outpri[1]};
+