From 0d361fd0a6c81a556d7b3b7fcee990ec451efde1 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Sat, 6 Mar 2010 14:22:00 +0000
Subject: [PATCH]
---
Common/LuaBindings.cpp | 2 +
Geo/GFaceCompound.cpp | 52 +--
Geo/MLine.cpp | 10 +
Geo/MLine.h | 1 +
Mesh/multiscalePartition.cpp | 85 ++--
Solver/CMakeLists.txt | 1 +
Solver/TESTCASES/Advection1D.lua | 13 +-
Solver/TESTCASES/ForwardFacingStep.lua | 2 +-
Solver/TESTCASES/ForwardFacingStepNEW.lua | 82 ++++
Solver/TESTCASES/MultirateAdvection.lua | 2 +-
Solver/TESTCASES/edge.geo | 2 +
Solver/TESTCASES/rect.geo | 2 +
Solver/TESTCASES/validation/Advection1D.lua | 75 ++++
Solver/TESTCASES/validation/edge.msh | 416 ++++++++++++++++++++
Solver/dgConservationLawShallowWater2d.cpp | 2 +-
Solver/dgConservationLawShallowWater2d.h | 2 +-
Solver/dgResidual.cpp | 7 +-
Solver/dgRungeKutta.cpp | 25 +-
Solver/dgRungeKutta.h | 2 +
Solver/eigenSolver.cpp | 2 +-
Solver/multiscaleLaplace.cpp | 21 +-
21 files changed, 719 insertions(+), 87 deletions(-)
create mode 100644 Solver/TESTCASES/ForwardFacingStepNEW.lua
create mode 100644 Solver/TESTCASES/validation/Advection1D.lua
create mode 100644 Solver/TESTCASES/validation/edge.msh
diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 8eba43f94c..b4f29dab2e 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -15,6 +15,7 @@
#include "dgGroupOfElements.h"
#include "dgDofContainer.h"
#include "dgConservationLawShallowWater2d.h"
+#include "dgConservationLawShallowWater1d.h"
#include "dgConservationLawAdvection.h"
#include "dgConservationLawPerfectGas.h"
#include "dgConservationLawWaveEquation.h"
@@ -344,6 +345,7 @@ binding::binding(){
dgConservationLaw::registerBindings(this);
dgConservationLawAdvectionDiffusionRegisterBindings(this);
dgConservationLawMaxwellRegisterBindings(this);
+ dgConservationLawShallowWater1dRegisterBindings(this);
dgConservationLawShallowWater2dRegisterBindings(this);
dgConservationLawWaveEquationRegisterBindings(this);
dgDofContainer::registerBindings(this);
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index a81c8ef47a..96bf32dde1 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -503,15 +503,14 @@ bool GFaceCompound::parametrize() const
else if (_mapping == CONFORMAL){
Msg::Debug("Parametrizing surface %d with 'conformal map'", tag());
fillNeumannBCS();
- bool withoutFolding = parametrize_conformal_spectral() ;
- printStuff();
- //exit(1);
+ //bool withoutFolding = parametrize_conformal_spectral() ;
+ bool withoutFolding = parametrize_conformal();
if ( withoutFolding == false ){
+ //printStuff(); exit(1);
Msg::Warning("$$$ Parametrization switched to harmonic map");
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
}
- //exit(1);
}
//Distance function
//-----------------
@@ -1138,7 +1137,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
cross21.addToMatrix(myAssembler, &se);
}
}
- double epsilon = 1.e-6;
+ double epsilon = 1.e-7;
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
if (std::find(ordered.begin(), ordered.end(), v) == ordered.end() ){
@@ -1167,14 +1166,19 @@ bool GFaceCompound::parametrize_conformal_spectral() const
myAssembler.assemble(v, 0, 2, v, 0, 2, 1.0);
}
}
+
// int NB = ordered.size();
-// for(std::vector<MVertex *>::iterator itv1 = ordered.begin(); itv1 !=ordered.end() ; ++itv1){
-// for(std::vector<MVertex *>::iterator itv2 = ordered.begin(); itv2 !=ordered.end() ; ++itv2){
-// myAssembler.assemble(*itv1, 0, 1, *itv2, 0, 1, -1/NB);
-// myAssembler.assemble(*itv1, 0, 2, *itv2, 0, 2, -1/NB);
+// for (int i = 0; i< NB; i++){
+// MVertex *v1 = ordered[i];
+// for (int j = i; j< NB; j++){
+// MVertex *v2 = ordered[j];
+// myAssembler.assemble(v1, 0, 1, v2, 0, 1, -1./NB);
+// myAssembler.assemble(v1, 0, 2, v2, 0, 2, -1./NB);
+// myAssembler.assemble(v2, 0, 1, v1, 0, 1, -1./NB);
+// myAssembler.assemble(v2, 0, 2, v1, 0, 2, -1./NB);
// }
// }
-
+
// diagBCTerm diag1(0, 1, &ONE);
// diagBCTerm diag2(0, 2, &ONE);
// it = _compound.begin();
@@ -1188,7 +1192,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
//-------------------------------
eigenSolver eig(&myAssembler, "B" , "A", true);
- bool converged = eig.solve(2, "largest");
+ bool converged = eig.solve(1, "largest");
if(converged) {
int k = 0;
@@ -1202,19 +1206,19 @@ bool GFaceCompound::parametrize_conformal_spectral() const
}
//if folding take second sallest eigenvalue
- bool noFolding = checkFolding(ordered);
- if (!noFolding ){
- coordinates.clear();
- int k = 0;
- std::vector<std::complex<double> > &ev = eig.getEigenVector(1);
- for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
- MVertex *v = *itv;
- double paramu = ev[k].real();
- double paramv = ev[k+1].real();
- coordinates[v] = SPoint3(paramu,paramv,0.0);
- k = k+2;
- }
- }
+// bool noFolding = checkFolding(ordered);
+// if (!noFolding ){
+// coordinates.clear();
+// int k = 0;
+// std::vector<std::complex<double> > &ev = eig.getEigenVector(1);
+// for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+// MVertex *v = *itv;
+// double paramu = ev[k].real();
+// double paramv = ev[k+1].real();
+// coordinates[v] = SPoint3(paramu,paramv,0.0);
+// k = k+2;
+// }
+// }
lsysA->clear();
lsysB->clear();
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 78bc74d741..77be232eb3 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -39,3 +39,13 @@ void MLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
*npts = nbP;
*pts = GQL;
}
+
+double MLine::getInnerRadius()
+{
+#if defined(HAVE_MESH)
+ double dist = _v[0]->distance(_v[1]);
+ return dist;
+#else
+ return 0.;
+#endif
+}
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 43ae3809f1..9b81989581 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -37,6 +37,7 @@ class MLine : public MElement {
virtual int getDim() const { return 1; }
virtual int getNumVertices() const { return 2; }
virtual MVertex *getVertex(int num){ return _v[num]; }
+ virtual double getInnerRadius(); // length of segment line
virtual void getVertexInfo(const MVertex * vertex, int &ithVertex) const
{
ithVertex = _v[0] == vertex ? 0 : 1;
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index 0fba0a7960..6d17375315 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -167,51 +167,52 @@ static int getAspectRatio(std::vector<MElement *> &elements,
AR = (int) ceil(2*3.14*area3D/(tot_length*tot_length));
}
-// std::set<MVertex*> vs;
-// for(unsigned int i = 0; i < elements.size(); i++){
-// MElement *e = elements[i];
-// for(int j = 0; j < e->getNumVertices(); j++){
-// vs.insert(e->getVertex(j));
-// }
-// }
-// SBoundingBox3d bb;
-// std::vector<SPoint3> vertices;
-// for (std::set<MVertex* >::iterator it = vs.begin(); it != vs.end(); it++){
-// SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
-// vertices.push_back(pt);
-// bb += pt;
-// }
-// double H = norm(SVector3(bb.max(), bb.min()));
-
-// //SOrientedBoundingBox obbox = SOrientedBoundingBox::buildOBB(vertices);
-// //double H = obbox.getMaxSize();
+ //compute AR also with Bounding box
+ std::set<MVertex*> vs;
+ for(unsigned int i = 0; i < elements.size(); i++){
+ MElement *e = elements[i];
+ for(int j = 0; j < e->getNumVertices(); j++){
+ vs.insert(e->getVertex(j));
+ }
+ }
+ SBoundingBox3d bb;
+ std::vector<SPoint3> vertices;
+ for (std::set<MVertex* >::iterator it = vs.begin(); it != vs.end(); it++){
+ SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
+ vertices.push_back(pt);
+ bb += pt;
+ }
+ double H = norm(SVector3(bb.max(), bb.min()));
+
+ //SOrientedBoundingBox obbox = SOrientedBoundingBox::buildOBB(vertices);
+ //double H = obbox.getMaxSize();
-// double D = H;
-// if (boundaries.size() > 0 ) D = 0.0;
-// for (unsigned int i = 0; i < boundaries.size(); i++){
-// std::set<MVertex*> vb;
-// std::vector<MEdge> iBound = boundaries[i];
-// for (unsigned int j = 0; j < iBound.size(); j++){
-// MEdge e = iBound[j];
-// vb.insert(e.getVertex(0));
-// vb.insert(e.getVertex(1));
-// }
-// std::vector<SPoint3> vBounds;
-// SBoundingBox3d bb;
-// for (std::set<MVertex* >::iterator it = vb.begin(); it != vb.end(); it++){
-// SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
-// vBounds.push_back(pt);
-// bb +=pt;
-// }
-// double iD = norm(SVector3(bb.max(), bb.min()));
-// D = std::max(D, iD);
+ double D = H;
+ if (boundaries.size() > 0 ) D = 0.0;
+ for (unsigned int i = 0; i < boundaries.size(); i++){
+ std::set<MVertex*> vb;
+ std::vector<MEdge> iBound = boundaries[i];
+ for (unsigned int j = 0; j < iBound.size(); j++){
+ MEdge e = iBound[j];
+ vb.insert(e.getVertex(0));
+ vb.insert(e.getVertex(1));
+ }
+ std::vector<SPoint3> vBounds;
+ SBoundingBox3d bb;
+ for (std::set<MVertex* >::iterator it = vb.begin(); it != vb.end(); it++){
+ SPoint3 pt((*it)->x(),(*it)->y(), (*it)->z());
+ vBounds.push_back(pt);
+ bb +=pt;
+ }
+ double iD = norm(SVector3(bb.max(), bb.min()));
+ D = std::max(D, iD);
-// //SOrientedBoundingBox obboxD = SOrientedBoundingBox::buildOBB(vBounds);
-// //D = std::max(D, obboxD.getMaxSize());
-// }
-// int AR = (int)ceil(H/D);
+ //SOrientedBoundingBox obboxD = SOrientedBoundingBox::buildOBB(vBounds);
+ //D = std::max(D, obboxD.getMaxSize());
+ }
+ int AR2 = (int)ceil(H/D);
- return AR;
+ return std::min(AR, AR2);
}
static void getGenusAndRatio(std::vector<MElement *> &elements, int & genus, int &AR, int &NB)
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index f5c4fc47a2..e88ccd212e 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -18,6 +18,7 @@ set(SRC
dgConservationLaw.cpp
dgConservationLawAdvection.cpp
dgSystemOfEquations.cpp
+ dgConservationLawShallowWater1d.cpp
dgConservationLawShallowWater2d.cpp
dgConservationLawWaveEquation.cpp
dgConservationLawMaxwell.cpp
diff --git a/Solver/TESTCASES/Advection1D.lua b/Solver/TESTCASES/Advection1D.lua
index d45996be69..c6f26566ef 100644
--- a/Solver/TESTCASES/Advection1D.lua
+++ b/Solver/TESTCASES/Advection1D.lua
@@ -63,12 +63,13 @@ solution:exportMsh('output/init_limit')
print'*** solve ***'
local x = os.clock()
-n = 5
dt = 0.03
-for i=1,100*n do
+for i=1,100 do
norm = rk:iterate44(law,dt,solution)
- if (i % n == 0) then
- print('|ITER|',i,'|NORM|',norm,'|DT|',dt,'|CPU|',os.clock() - x)
- solution:exportMsh(string.format('output/solution-%06d',i))
- end
+ if (i % 20 == 0) then
+ print('|ITER|',i,'|NORM|',norm,'|DT|',dt,'|CPU|',os.clock() - x)
+ end
+ if (i % 200 == 0) then
+ solution:exportMsh(string.format('output/solution-%06d',i))
+ end
end
diff --git a/Solver/TESTCASES/ForwardFacingStep.lua b/Solver/TESTCASES/ForwardFacingStep.lua
index 0507416ca7..ebe99da4e9 100644
--- a/Solver/TESTCASES/ForwardFacingStep.lua
+++ b/Solver/TESTCASES/ForwardFacingStep.lua
@@ -41,7 +41,7 @@ FS = functionLua(4, 'free_stream', {'XYZ'}):getName()
law=dgPerfectGasLaw2d()
DG:setConservationLaw(law)
-law:addBoundaryCondition('Walls',law:newBoundaryWall())
+law:addBoundaryCondition('Walls',law:newNonSlipWallBoundary())
law:addBoundaryCondition('LeftRight',law:newOutsideValueBoundary(FS))
--law:addBoundaryCondition('Walls',law:newOutsideValueBoundary(FS))
diff --git a/Solver/TESTCASES/ForwardFacingStepNEW.lua b/Solver/TESTCASES/ForwardFacingStepNEW.lua
new file mode 100644
index 0000000000..1eb870f462
--- /dev/null
+++ b/Solver/TESTCASES/ForwardFacingStepNEW.lua
@@ -0,0 +1,82 @@
+MACH = 3.0;
+GAMMA = 1.4;
+U = 3.0
+V = 0.0
+RHO = 1.4;
+
+PRES = RHO*U*U/(GAMMA*MACH*MACH)
+
+--PRES = 1;
+--PRES = ./(MACH*RHO*RHO*GAMMA*GAMMA)
+
+SOUND = math.sqrt(U*U+V*V)/MACH
+
+--[[
+ Function for initial conditions
+--]]
+function free_stream( XYZ, FCT )
+ for i=0,XYZ:size1()-1 do
+ FCT:set(i,0,RHO)
+ FCT:set(i,1,RHO*U)
+ FCT:set(i,2,RHO*V)
+ FCT:set(i,3, 0.5*RHO*(U*U+V*V)+PRES/(GAMMA-1))
+ end
+end
+
+--[[
+ Example of a lua program driving the DG code
+--]]
+
+
+print'*** Loading the mesh and the model ***'
+model = GModel ()
+model:load ('step.geo')
+model:load ('step.msh')
+order = 1
+dimension = 2
+
+FS = functionLua(4, 'free_stream', {'XYZ'}):getName()
+
+-- boundary condition
+law=dgPerfectGasLaw2d()
+law:addBoundaryCondition('Walls',law:newNonSlipWallBoundary())
+law:addBoundaryCondition('LeftRight',law:newOutsideValueBoundary(FS))
+--law:addBoundaryCondition('Walls',law:newOutsideValueBoundary(FS))
+
+groups = dgGroupCollection(model, dimension, order)
+groups:buildGroupsOfInterfaces()
+
+-- build Runge Kutta and limiter
+rk=dgRungeKutta()
+limiter = dgSlopeLimiter(law)
+rk:setLimiter(limiter)
+
+-- build solution vector
+solution = dgDofContainer(groups, law:getNbFields())
+solution:L2Projection(FS)
+solution:exportGroupIdMsh()
+
+print'*** setting the initial solution ***'
+solution:exportMsh('output/init')
+limiter:apply(solution)
+solution:exportMsh('output/init_limit')
+
+print'*** solve ***'
+dg = dgSystemOfEquations (model)
+CFL = 2
+local x = os.clock()
+for i=1,5000 do
+ dt = CFL * rk:computeInvSpectralRadius(law,solution);
+-- norm = rk:iterate44(law,dt,solution)
+ norm = rk:iterateEuler(law,dt,solution)
+ if (i % 10 == 0) then
+ print('|ITER|',i,'|NORM|',norm,'|DT|',dt,'|CPU|',os.clock() - x)
+ end
+ if (i % 100 == 0) then
+ solution:exportMsh(string.format('output/solution-%06d', i))
+ end
+end
+
+print'*** done ***'
+
+
diff --git a/Solver/TESTCASES/MultirateAdvection.lua b/Solver/TESTCASES/MultirateAdvection.lua
index 47636aa8ba..7e1e0a4184 100644
--- a/Solver/TESTCASES/MultirateAdvection.lua
+++ b/Solver/TESTCASES/MultirateAdvection.lua
@@ -60,7 +60,7 @@ FS = functionLua(1, 'initial_condition', {'XYZ'}):getName()
GC=dgGroupCollection(myModel,2,order)
solTmp=dgDofContainer(GC,1)
solTmp:L2Projection(FS)
-dt=GC:splitGroupsForMultirate(100,law,solTmp)
+dt=GC:splitGroupsForMultirate(4,law,solTmp)
GC:buildGroupsOfInterfaces(myModel,2,order)
solution=dgDofContainer(GC,1)
solution:L2Projection(FS)
diff --git a/Solver/TESTCASES/edge.geo b/Solver/TESTCASES/edge.geo
index 67fd0dd760..6c17d4e77a 100644
--- a/Solver/TESTCASES/edge.geo
+++ b/Solver/TESTCASES/edge.geo
@@ -1,6 +1,8 @@
+
Point(1) = {-1, 0, 0, 0.01};
Point(2) = {1, 0, 0, 0.01};
Line(1) = {1, 2};
+
Physical Point("Left") = {1};
Physical Point("Right") = {2};
Physical Line("Line") = {1};
diff --git a/Solver/TESTCASES/rect.geo b/Solver/TESTCASES/rect.geo
index ae5fbc2344..c124edd203 100644
--- a/Solver/TESTCASES/rect.geo
+++ b/Solver/TESTCASES/rect.geo
@@ -15,9 +15,11 @@ Line Loop(8) = {7, -2, 5, 6};
Plane Surface(9) = {8};
Line Loop(10) = {3, 4, 1, 2};
Plane Surface(11) = {10};
+
Physical Surface("sprut") = {11, 9};
Physical Line("Walls") = {5, 7, 3, 4, 1};
Physical Line("Top") = {6};
+
Recombine Surface {9, 11};
Field[1] = MathEval;
diff --git a/Solver/TESTCASES/validation/Advection1D.lua b/Solver/TESTCASES/validation/Advection1D.lua
new file mode 100644
index 0000000000..c6f26566ef
--- /dev/null
+++ b/Solver/TESTCASES/validation/Advection1D.lua
@@ -0,0 +1,75 @@
+
+--[[
+ Function for initial conditions
+--]]
+
+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)
+ if (x>-0.3 and x<0.3) then
+ f:set (i, 0, 1)
+ else
+ f:set (i, 0, 0)
+ end
+ end
+end
+
+--[[
+ Example of a lua program driving the DG code
+--]]
+
+model = GModel ()
+model:load ('edge.msh')
+order=1
+dimension=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)
+
+law = dgConservationLawAdvectionDiffusion(functionConstant(v):getName(), '')
+
+-- boundary condition
+outside=fullMatrix(1,1)
+outside:set(0,0,0.)
+bndcondition=law:newOutsideValueBoundary(functionConstant(outside):getName())
+law:addBoundaryCondition('Left',bndcondition)
+law:addBoundaryCondition('Right',bndcondition)
+
+groups = dgGroupCollection(model, dimension, order)
+groups:buildGroupsOfInterfaces()
+
+-- build Runge Kutta and limiter
+rk=dgRungeKutta()
+limiter = dgSlopeLimiter(law)
+rk:setLimiter(limiter)
+
+-- build solution vector
+FS = functionLua(1, 'initial_condition', {'XYZ'}):getName()
+solution = dgDofContainer(groups, law:getNbFields())
+solution:L2Projection(FS)
+
+solution:exportMsh('output/init')
+limiter:apply(solution)
+solution:exportMsh('output/init_limit')
+
+print'*** solve ***'
+local x = os.clock()
+dt = 0.03
+for i=1,100 do
+ norm = rk:iterate44(law,dt,solution)
+ if (i % 20 == 0) then
+ print('|ITER|',i,'|NORM|',norm,'|DT|',dt,'|CPU|',os.clock() - x)
+ end
+ if (i % 200 == 0) then
+ solution:exportMsh(string.format('output/solution-%06d',i))
+ end
+end
diff --git a/Solver/TESTCASES/validation/edge.msh b/Solver/TESTCASES/validation/edge.msh
new file mode 100644
index 0000000000..b77e58a920
--- /dev/null
+++ b/Solver/TESTCASES/validation/edge.msh
@@ -0,0 +1,416 @@
+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$PhysicalNames
+3
+0 1 "Left"
+0 2 "Right"
+1 3 "Line"
+$EndPhysicalNames
+$Nodes
+200
+1 -1 0 0
+2 1 0 0
+3 -0.9899497487437465 0 0
+4 -0.9798994974874929 0 0
+5 -0.9698492462312395 0 0
+6 -0.959798994974986 0 0
+7 -0.9497487437187324 0 0
+8 -0.9396984924624789 0 0
+9 -0.9296482412062255 0 0
+10 -0.9195979899499719 0 0
+11 -0.9095477386937184 0 0
+12 -0.8994974874374648 0 0
+13 -0.8894472361812114 0 0
+14 -0.8793969849249579 0 0
+15 -0.8693467336687044 0 0
+16 -0.8592964824124508 0 0
+17 -0.8492462311561974 0 0
+18 -0.8391959798999438 0 0
+19 -0.8291457286436903 0 0
+20 -0.8190954773874368 0 0
+21 -0.8090452261311832 0 0
+22 -0.7989949748749298 0 0
+23 -0.7889447236186762 0 0
+24 -0.7788944723624227 0 0
+25 -0.7688442211061692 0 0
+26 -0.7587939698499158 0 0
+27 -0.7487437185936622 0 0
+28 -0.7386934673374087 0 0
+29 -0.7286432160811551 0 0
+30 -0.7185929648249016 0 0
+31 -0.7085427135686482 0 0
+32 -0.6984924623123947 0 0
+33 -0.6884422110561411 0 0
+34 -0.6783919597998875 0 0
+35 -0.6683417085436341 0 0
+36 -0.6582914572873806 0 0
+37 -0.6482412060311271 0 0
+38 -0.6381909547748736 0 0
+39 -0.6281407035186201 0 0
+40 -0.6180904522623665 0 0
+41 -0.608040201006113 0 0
+42 -0.5979899497498595 0 0
+43 -0.5879396984936061 0 0
+44 -0.5778894472373525 0 0
+45 -0.567839195981099 0 0
+46 -0.5577889447248454 0 0
+47 -0.5477386934685919 0 0
+48 -0.5376884422123385 0 0
+49 -0.527638190956085 0 0
+50 -0.5175879396998314 0 0
+51 -0.5075376884435778 0 0
+52 -0.4974874371873209 0 0
+53 -0.4874371859310535 0 0
+54 -0.477386934674786 0 0
+55 -0.4673366834185185 0 0
+56 -0.4572864321622511 0 0
+57 -0.4472361809059836 0 0
+58 -0.4371859296497161 0 0
+59 -0.4271356783934487 0 0
+60 -0.4170854271371812 0 0
+61 -0.4070351758809138 0 0
+62 -0.3969849246246463 0 0
+63 -0.3869346733683788 0 0
+64 -0.3768844221121114 0 0
+65 -0.3668341708558439 0 0
+66 -0.3567839195995764 0 0
+67 -0.346733668343309 0 0
+68 -0.3366834170870415 0 0
+69 -0.3266331658307741 0 0
+70 -0.3165829145745066 0 0
+71 -0.3065326633182391 0 0
+72 -0.2964824120619717 0 0
+73 -0.2864321608057042 0 0
+74 -0.2763819095494368 0 0
+75 -0.2663316582931693 0 0
+76 -0.2562814070369018 0 0
+77 -0.2462311557806344 0 0
+78 -0.2361809045243669 0 0
+79 -0.2261306532680994 0 0
+80 -0.216080402011832 0 0
+81 -0.2060301507555645 0 0
+82 -0.1959798994992971 0 0
+83 -0.1859296482430296 0 0
+84 -0.1758793969867621 0 0
+85 -0.1658291457304947 0 0
+86 -0.1557788944742272 0 0
+87 -0.1457286432179598 0 0
+88 -0.1356783919616923 0 0
+89 -0.1256281407054248 0 0
+90 -0.1155778894491574 0 0
+91 -0.1055276381928899 0 0
+92 -0.09547738693662244 0 0
+93 -0.08542713568035498 0 0
+94 -0.07537688442408752 0 0
+95 -0.06532663316782006 0 0
+96 -0.0552763819115526 0 0
+97 -0.04522613065528525 0 0
+98 -0.03517587939901778 0 0
+99 -0.02512562814275032 0 0
+100 -0.01507537688648286 0 0
+101 -0.005025125630215399 0 0
+102 0.005025125626066052 0 0
+103 0.01507537688236149 0 0
+104 0.02512562813865671 0 0
+105 0.03517587939495215 0 0
+106 0.04522613065124759 0 0
+107 0.0552763819075428 0 0
+108 0.06532663316383824 0 0
+109 0.07537688442013346 0 0
+110 0.0854271356764289 0 0
+111 0.09547738693272434 0 0
+112 0.1055276381890196 0 0
+113 0.115577889445315 0 0
+114 0.1256281407016102 0 0
+115 0.1356783919579057 0 0
+116 0.1457286432142011 0 0
+117 0.1557788944704963 0 0
+118 0.1658291457267917 0 0
+119 0.1758793969830872 0 0
+120 0.1859296482393824 0 0
+121 0.1959798994956778 0 0
+122 0.2060301507519731 0 0
+123 0.2160804020082685 0 0
+124 0.2261306532645639 0 0
+125 0.2361809045208592 0 0
+126 0.2462311557771546 0 0
+127 0.25628140703345 0 0
+128 0.2663316582897453 0 0
+129 0.2763819095460407 0 0
+130 0.2864321608023359 0 0
+131 0.2964824120586314 0 0
+132 0.3065326633149268 0 0
+133 0.316582914571222 0 0
+134 0.3266331658275174 0 0
+135 0.3366834170838127 0 0
+136 0.3467336683401081 0 0
+137 0.3567839195964035 0 0
+138 0.3668341708526988 0 0
+139 0.3768844221089942 0 0
+140 0.3869346733652896 0 0
+141 0.3969849246215849 0 0
+142 0.4070351758778803 0 0
+143 0.4170854271341755 0 0
+144 0.427135678390471 0 0
+145 0.4371859296467664 0 0
+146 0.4472361809030616 0 0
+147 0.457286432159357 0 0
+148 0.4673366834156525 0 0
+149 0.4773869346719477 0 0
+150 0.4874371859282431 0 0
+151 0.4974874371845384 0 0
+152 0.5075376884408442 0 0
+153 0.5175879396971537 0 0
+154 0.5276381909534629 0 0
+155 0.5376884422097721 0 0
+156 0.5477386934660815 0 0
+157 0.5577889447223907 0 0
+158 0.5678391959787001 0 0
+159 0.5778894472350093 0 0
+160 0.5879396984913188 0 0
+161 0.597989949747628 0 0
+162 0.6080402010039372 0 0
+163 0.6180904522602466 0 0
+164 0.6281407035165558 0 0
+165 0.6381909547728652 0 0
+166 0.6482412060291745 0 0
+167 0.6582914572854837 0 0
+168 0.6683417085417931 0 0
+169 0.6783919597981023 0 0
+170 0.6884422110544117 0 0
+171 0.6984924623107209 0 0
+172 0.7085427135670304 0 0
+173 0.7185929648233396 0 0
+174 0.7286432160796488 0 0
+175 0.7386934673359582 0 0
+176 0.7487437185922674 0 0
+177 0.7587939698485768 0 0
+178 0.768844221104886 0 0
+179 0.7788944723611955 0 0
+180 0.7889447236175047 0 0
+181 0.7989949748738139 0 0
+182 0.8090452261301233 0 0
+183 0.8190954773864325 0 0
+184 0.8291457286427419 0 0
+185 0.8391959798990511 0 0
+186 0.8492462311553606 0 0
+187 0.8592964824116698 0 0
+188 0.869346733667979 0 0
+189 0.8793969849242884 0 0
+190 0.8894472361805976 0 0
+191 0.8994974874369071 0 0
+192 0.9095477386932163 0 0
+193 0.9195979899495255 0 0
+194 0.9296482412058349 0 0
+195 0.9396984924621441 0 0
+196 0.9497487437184535 0 0
+197 0.9597989949747627 0 0
+198 0.9698492462310722 0 0
+199 0.9798994974873814 0 0
+200 0.9899497487436906 0 0
+$EndNodes
+$Elements
+201
+1 15 2 1 1 1
+2 15 2 2 2 2
+3 1 2 3 1 1 3
+4 1 2 3 1 3 4
+5 1 2 3 1 4 5
+6 1 2 3 1 5 6
+7 1 2 3 1 6 7
+8 1 2 3 1 7 8
+9 1 2 3 1 8 9
+10 1 2 3 1 9 10
+11 1 2 3 1 10 11
+12 1 2 3 1 11 12
+13 1 2 3 1 12 13
+14 1 2 3 1 13 14
+15 1 2 3 1 14 15
+16 1 2 3 1 15 16
+17 1 2 3 1 16 17
+18 1 2 3 1 17 18
+19 1 2 3 1 18 19
+20 1 2 3 1 19 20
+21 1 2 3 1 20 21
+22 1 2 3 1 21 22
+23 1 2 3 1 22 23
+24 1 2 3 1 23 24
+25 1 2 3 1 24 25
+26 1 2 3 1 25 26
+27 1 2 3 1 26 27
+28 1 2 3 1 27 28
+29 1 2 3 1 28 29
+30 1 2 3 1 29 30
+31 1 2 3 1 30 31
+32 1 2 3 1 31 32
+33 1 2 3 1 32 33
+34 1 2 3 1 33 34
+35 1 2 3 1 34 35
+36 1 2 3 1 35 36
+37 1 2 3 1 36 37
+38 1 2 3 1 37 38
+39 1 2 3 1 38 39
+40 1 2 3 1 39 40
+41 1 2 3 1 40 41
+42 1 2 3 1 41 42
+43 1 2 3 1 42 43
+44 1 2 3 1 43 44
+45 1 2 3 1 44 45
+46 1 2 3 1 45 46
+47 1 2 3 1 46 47
+48 1 2 3 1 47 48
+49 1 2 3 1 48 49
+50 1 2 3 1 49 50
+51 1 2 3 1 50 51
+52 1 2 3 1 51 52
+53 1 2 3 1 52 53
+54 1 2 3 1 53 54
+55 1 2 3 1 54 55
+56 1 2 3 1 55 56
+57 1 2 3 1 56 57
+58 1 2 3 1 57 58
+59 1 2 3 1 58 59
+60 1 2 3 1 59 60
+61 1 2 3 1 60 61
+62 1 2 3 1 61 62
+63 1 2 3 1 62 63
+64 1 2 3 1 63 64
+65 1 2 3 1 64 65
+66 1 2 3 1 65 66
+67 1 2 3 1 66 67
+68 1 2 3 1 67 68
+69 1 2 3 1 68 69
+70 1 2 3 1 69 70
+71 1 2 3 1 70 71
+72 1 2 3 1 71 72
+73 1 2 3 1 72 73
+74 1 2 3 1 73 74
+75 1 2 3 1 74 75
+76 1 2 3 1 75 76
+77 1 2 3 1 76 77
+78 1 2 3 1 77 78
+79 1 2 3 1 78 79
+80 1 2 3 1 79 80
+81 1 2 3 1 80 81
+82 1 2 3 1 81 82
+83 1 2 3 1 82 83
+84 1 2 3 1 83 84
+85 1 2 3 1 84 85
+86 1 2 3 1 85 86
+87 1 2 3 1 86 87
+88 1 2 3 1 87 88
+89 1 2 3 1 88 89
+90 1 2 3 1 89 90
+91 1 2 3 1 90 91
+92 1 2 3 1 91 92
+93 1 2 3 1 92 93
+94 1 2 3 1 93 94
+95 1 2 3 1 94 95
+96 1 2 3 1 95 96
+97 1 2 3 1 96 97
+98 1 2 3 1 97 98
+99 1 2 3 1 98 99
+100 1 2 3 1 99 100
+101 1 2 3 1 100 101
+102 1 2 3 1 101 102
+103 1 2 3 1 102 103
+104 1 2 3 1 103 104
+105 1 2 3 1 104 105
+106 1 2 3 1 105 106
+107 1 2 3 1 106 107
+108 1 2 3 1 107 108
+109 1 2 3 1 108 109
+110 1 2 3 1 109 110
+111 1 2 3 1 110 111
+112 1 2 3 1 111 112
+113 1 2 3 1 112 113
+114 1 2 3 1 113 114
+115 1 2 3 1 114 115
+116 1 2 3 1 115 116
+117 1 2 3 1 116 117
+118 1 2 3 1 117 118
+119 1 2 3 1 118 119
+120 1 2 3 1 119 120
+121 1 2 3 1 120 121
+122 1 2 3 1 121 122
+123 1 2 3 1 122 123
+124 1 2 3 1 123 124
+125 1 2 3 1 124 125
+126 1 2 3 1 125 126
+127 1 2 3 1 126 127
+128 1 2 3 1 127 128
+129 1 2 3 1 128 129
+130 1 2 3 1 129 130
+131 1 2 3 1 130 131
+132 1 2 3 1 131 132
+133 1 2 3 1 132 133
+134 1 2 3 1 133 134
+135 1 2 3 1 134 135
+136 1 2 3 1 135 136
+137 1 2 3 1 136 137
+138 1 2 3 1 137 138
+139 1 2 3 1 138 139
+140 1 2 3 1 139 140
+141 1 2 3 1 140 141
+142 1 2 3 1 141 142
+143 1 2 3 1 142 143
+144 1 2 3 1 143 144
+145 1 2 3 1 144 145
+146 1 2 3 1 145 146
+147 1 2 3 1 146 147
+148 1 2 3 1 147 148
+149 1 2 3 1 148 149
+150 1 2 3 1 149 150
+151 1 2 3 1 150 151
+152 1 2 3 1 151 152
+153 1 2 3 1 152 153
+154 1 2 3 1 153 154
+155 1 2 3 1 154 155
+156 1 2 3 1 155 156
+157 1 2 3 1 156 157
+158 1 2 3 1 157 158
+159 1 2 3 1 158 159
+160 1 2 3 1 159 160
+161 1 2 3 1 160 161
+162 1 2 3 1 161 162
+163 1 2 3 1 162 163
+164 1 2 3 1 163 164
+165 1 2 3 1 164 165
+166 1 2 3 1 165 166
+167 1 2 3 1 166 167
+168 1 2 3 1 167 168
+169 1 2 3 1 168 169
+170 1 2 3 1 169 170
+171 1 2 3 1 170 171
+172 1 2 3 1 171 172
+173 1 2 3 1 172 173
+174 1 2 3 1 173 174
+175 1 2 3 1 174 175
+176 1 2 3 1 175 176
+177 1 2 3 1 176 177
+178 1 2 3 1 177 178
+179 1 2 3 1 178 179
+180 1 2 3 1 179 180
+181 1 2 3 1 180 181
+182 1 2 3 1 181 182
+183 1 2 3 1 182 183
+184 1 2 3 1 183 184
+185 1 2 3 1 184 185
+186 1 2 3 1 185 186
+187 1 2 3 1 186 187
+188 1 2 3 1 187 188
+189 1 2 3 1 188 189
+190 1 2 3 1 189 190
+191 1 2 3 1 190 191
+192 1 2 3 1 191 192
+193 1 2 3 1 192 193
+194 1 2 3 1 193 194
+195 1 2 3 1 194 195
+196 1 2 3 1 195 196
+197 1 2 3 1 196 197
+198 1 2 3 1 197 198
+199 1 2 3 1 198 199
+200 1 2 3 1 199 200
+201 1 2 3 1 200 2
+$EndElements
diff --git a/Solver/dgConservationLawShallowWater2d.cpp b/Solver/dgConservationLawShallowWater2d.cpp
index 572c219ac2..6f742b4a75 100644
--- a/Solver/dgConservationLawShallowWater2d.cpp
+++ b/Solver/dgConservationLawShallowWater2d.cpp
@@ -176,7 +176,7 @@ dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryWall(){
#include "Bindings.h"
void dgConservationLawShallowWater2dRegisterBindings (binding *b){
classBinding *cb = b->addClass<dgConservationLawShallowWater2d>("dgConservationLawShallowWater2d");
- cb->setDescription("The conservative sallow water conservation law. (H, Hu, Hv)");
+ cb->setDescription("The non-conservative shallow water conservation law. (eta, u, v)");
methodBinding *cm;
cm = cb->addMethod("newBoundaryWall",&dgConservationLawShallowWater2d::newBoundaryWall);
cm->setDescription("slip wall boundary");
diff --git a/Solver/dgConservationLawShallowWater2d.h b/Solver/dgConservationLawShallowWater2d.h
index 2cf268972b..404944542d 100644
--- a/Solver/dgConservationLawShallowWater2d.h
+++ b/Solver/dgConservationLawShallowWater2d.h
@@ -17,7 +17,7 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const;
dgConservationLawShallowWater2d()
{
- _nbf = 3; // H U(=Hu) V(=Hv)
+ _nbf = 3; // eta u v
}
dgBoundaryCondition *newBoundaryWall();
};
diff --git a/Solver/dgResidual.cpp b/Solver/dgResidual.cpp
index 0289ba87d3..57005f39e7 100644
--- a/Solver/dgResidual.cpp
+++ b/Solver/dgResidual.cpp
@@ -50,10 +50,10 @@ void dgResidualVolume::compute1Group(dgGroupOfElements &group, fullMatrix<double
fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * _nbFields)};
fullMatrix<double> fuvwe;
fullMatrix<double> source;
- fullMatrix<double> dPsiDx,dofs;
+ fullMatrix<double> dPsiDx,dofs;
// ----- 2.3 --- iterate on elements
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
- // ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
+ // ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
_solutionQPe.setAsProxy(solutionQP, iElement*_nbFields, _nbFields );
if(_gradientSolutionQPe.somethingDependOnMe()){
@@ -97,13 +97,12 @@ void dgResidualVolume::compute1Group(dgGroupOfElements &group, fullMatrix<double
}
}
}
-
// ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
if(_convectiveFlux || _diffusiveFlux){
for (int iUVW=0;iUVW<group.getDimUVW();iUVW++){
residual.gemm(group.getFluxRedistributionMatrix(iUVW),Fuvw[iUVW]);
}
- }
+ }
if(_sourceTerm){
residual.gemm(group.getSourceRedistributionMatrix(),Source);
}
diff --git a/Solver/dgRungeKutta.cpp b/Solver/dgRungeKutta.cpp
index 2c24bc36fd..63e07b2e34 100644
--- a/Solver/dgRungeKutta.cpp
+++ b/Solver/dgRungeKutta.cpp
@@ -1,11 +1,14 @@
#include <stdlib.h>
#include "dgRungeKutta.h"
+#include "dgAlgorithm.h"
#include "dgConservationLaw.h"
#include "dgDofContainer.h"
#include "dgLimiter.h"
#include "dgTransformNodalValue.h"
#include "dgResidual.h"
#include "dgGroupOfElements.h"
+#include <limits.h>
+#include <stdio.h>
double dgRungeKutta::iterateEuler(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
double A[] = {0};
@@ -211,6 +214,24 @@ double dgRungeKutta::nonDiagonalRK(const dgConservationLaw *claw,
return solution->norm();
}
+double dgRungeKutta::computeInvSpectralRadius(const dgConservationLaw *claw,
+ dgDofContainer *solution){
+ double sr = 1.e22;
+ dgGroupCollection *groups=solution->getGroups();
+ for (int i=0;i<groups->getNbElementGroups();i++){
+ std::vector<double> DTS;
+ dgAlgorithm::computeElementaryTimeSteps(*claw, *(groups->getElementGroup(i)), solution->getGroupProxy(i), DTS);
+ for (int k=0;k<DTS.size();k++) sr = std::min(sr,DTS[k]);
+ }
+ #ifdef HAVE_MPI
+ double sr_min;
+ MPI_Allreduce((void *)&sr, &sr_min, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+ return sr_min;
+ #else
+ return sr;
+ #endif
+}
+
dgRungeKutta::dgRungeKutta():_limiter(NULL),_TransformNodalValue(NULL){}
@@ -247,7 +268,9 @@ void dgRungeKutta::registerBindings(binding *b) {
cm = cb->addMethod("setLimiter",&dgRungeKutta::setLimiter);
cm->setArgNames("limiter",NULL);
cm->setDescription("if a limiter is set, it is applied after each RK stage");
-
+ cm = cb->addMethod("computeInvSpectralRadius",&dgRungeKutta::computeInvSpectralRadius);
+ cm->setArgNames("law","solution",NULL);
+ cm->setDescription("Returns the inverse of the spectral radius of L(U). Useful for computing stable explicit time step.");
cm = cb->addMethod("setTransformNodalValue",&dgRungeKutta::setTransformNodalValue);
cm->setArgNames("TransformNodalValue",NULL);
cm->setDescription("if the Nodal values is transformed in first step of RK");
diff --git a/Solver/dgRungeKutta.h b/Solver/dgRungeKutta.h
index 8088e783c8..4c178b7499 100644
--- a/Solver/dgRungeKutta.h
+++ b/Solver/dgRungeKutta.h
@@ -5,6 +5,7 @@
class dgConservationLaw;
class dgDofContainer;
class dgLimiter;
+class dgAlgorithm;
class dgTransformNodalValue;
class dgGroupCollection;
class dgGroupOfElements;
@@ -21,6 +22,7 @@ class dgRungeKutta {
dgTransformNodalValue *_TransformNodalValue;
public:
void setLimiter(dgLimiter *limiter) { _limiter = limiter; }
+ double computeInvSpectralRadius(const dgConservationLaw *claw, dgDofContainer *solution);
void setTransformNodalValue(dgTransformNodalValue *TransformNodalValue) { _TransformNodalValue = TransformNodalValue; }
double iterateEuler(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
diff --git a/Solver/eigenSolver.cpp b/Solver/eigenSolver.cpp
index 3568215d56..57f3534520 100644
--- a/Solver/eigenSolver.cpp
+++ b/Solver/eigenSolver.cpp
@@ -52,7 +52,7 @@ bool eigenSolver::solve(int numEigenValues, std::string which)
// set some default options
_try(EPSSetDimensions(eps, numEigenValues, PETSC_DECIDE, PETSC_DECIDE));
- _try(EPSSetTolerances(eps,1.e-7,100));
+ _try(EPSSetTolerances(eps,1.e-7,50));
_try(EPSSetType(eps,EPSARNOLDI)); //EPSKRYLOVSCHUR is default
//_try(EPSSetType(eps,EPSARPACK));
//_try(EPSSetType(eps,EPSPOWER));
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 2f8699976d..aef8830e52 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -227,6 +227,7 @@ static void recur_compute_centers_ (double R, double a1, double a2,
centers.push_back(std::make_pair(PL,zero));
centers.push_back(std::make_pair(PR,zero));
+ std::map<SPoint2, double> CR;
for (int i=0;i<root->children.size();i++){
multiscaleLaplaceLevel* m = root->children[i];
centers.push_back(std::make_pair(m->center,m));
@@ -237,12 +238,13 @@ static void recur_compute_centers_ (double R, double a1, double a2,
m->radius = std::max(m->radius,sqrt ((m->center.x() - p.x())*(m->center.x() - p.x())+
(m->center.y() - p.y())*(m->center.y() - p.y())));
}
+ CR.insert(std::make_pair(m->center, m->radius));
}
//add the center of real holes ...
- if (root->children.size()==0 ){//|| boundaries.size()-1 != root->children.size()){
- std::vector<std::vector<MEdge> > boundaries;
- connected_bounds(root->elements, boundaries);
+ std::vector<std::vector<MEdge> > boundaries;
+ connected_bounds(root->elements, boundaries);
+ if (root->children.size()==0 ){ // || boundaries.size()-1 != root->children.size()){
for (unsigned int i = 0; i < boundaries.size(); i++){
std::vector<MEdge> me = boundaries[i];
SPoint2 c(0.0,0.0);
@@ -258,6 +260,15 @@ static void recur_compute_centers_ (double R, double a1, double a2,
rad = std::max(rad,sqrt ((c.x() - p.x())*(c.x() - p.x())+
(c.y() - p.y())*(c.y() - p.y())));
}
+
+ //check if the center has not been added
+// bool newCenter = false;
+// for (std::map<SPoint2,double>::iterator it = CR.begin(); it != CR.end(); it++){
+// SPoint2 p = it->first;
+// double radius = it->second;
+// double dist = sqrt ((c.x() - p.x())*(c.x() - p.x())+ (c.y() - p.y())*(c.y() - p.y()));
+// if (dist > radius) newCenter = true;
+// }
if (std::abs(rad/root->radius) < 0.9 && std::abs(rad) < 0.99){
centers.push_back(std::make_pair(c,zero));
}
@@ -903,7 +914,7 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
MElement *e = level.elements[i];
std::vector<SPoint2> localCoord;
double local_size = localSize(e,solution);
- if (local_size < 1.e-5 * global_size)
+ if (local_size < 5.e-5 * global_size) //1.e-5
tooSmall.push_back(e);
else goodSize.push_back(e);
}
@@ -939,7 +950,7 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
for (int k=0; k<regions_[i].size() ; k++){
MElement *e = regions_[i][k];
double local_size = localSize(e,solution);
- if (local_size < 1.e-7 * global_size)
+ if (local_size < 1.e-7 * global_size) //1.e-7
really_small_elements = true;
}
if(really_small_elements ) regions.push_back(regions_[i]);
--
GitLab