From 088f5e6e558f9450e072ac1c5d76941132bee2f4 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Wed, 17 Mar 2010 16:47:28 +0000
Subject: [PATCH] Added bathymetry and dissipation to shallow 1d
---
Geo/GFaceCompound.cpp | 9 +-
Mesh/meshGFace.cpp | 6 +-
Mesh/multiscalePartition.cpp | 9 +-
Solver/TESTCASES/Aorta.lua | 67 ++++
Solver/TESTCASES/CylinderEddies.lua | 4 +-
Solver/TESTCASES/DamBreak.lua | 21 +-
Solver/TESTCASES/DamBreak2d.lua | 72 ++++
Solver/TESTCASES/ForwardFacingStepNEW.lua | 2 +-
Solver/TESTCASES/aorta.geo | 8 +
Solver/TESTCASES/edge.msh | 402 +++++++++++++++++++--
Solver/dgConservationLawShallowWater1d.cpp | 155 +++++---
Solver/dgConservationLawShallowWater1d.h | 23 +-
Solver/dgConservationLawShallowWater2d.cpp | 32 ++
Solver/multiscaleLaplace.cpp | 33 +-
benchmarks/3d/Torus.geo | 4 +-
15 files changed, 701 insertions(+), 146 deletions(-)
create mode 100644 Solver/TESTCASES/Aorta.lua
create mode 100644 Solver/TESTCASES/DamBreak2d.lua
create mode 100644 Solver/TESTCASES/aorta.geo
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index d3063533e8..012c4cd8a5 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -534,7 +534,7 @@ bool GFaceCompound::parametrize() const
double AR = checkAspectRatio();
if (floor(AR) > AR_MAX){
- Msg::Warning("Geometrical aspect ratio too high %d ", AR);
+ Msg::Warning("Geometrical aspect ratio too high AR=%d ", (int)AR);
//exit(1);
paramOK = true; //false;
}
@@ -1890,12 +1890,7 @@ double GFaceCompound::checkAspectRatio() const
double AR = 2*3.14*area3D/(tot_length*tot_length);
if (areaMin < limit && nb > 2) {
- Msg::Warning("Geometrical aspect ratio too high (a_2D=%g)", areaMin);
- SBoundingBox3d bboxH = bounds();
- double H = getSizeH();
- double D = getSizeBB(_U0);
- double eta = H/D;
- int nbSplit = -2;
+ Msg::Warning("Too small triangles in mapping (a_2D=%g)", areaMin);
}
else {
Msg::Debug("Geometrical aspect ratio is OK :-)");
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index cd518f885a..3b2971cbed 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1372,7 +1372,7 @@ void partitionAndRemesh(GFaceCompound *gf)
gf->model()->createTopologyFromFaces(pFaces);
Msg::Info("Multiscale Partition SUCCESSFULLY PERFORMED : %d parts", NF );
- gf->model()->writeMSH("multiscalePARTS.msh", 2.0, false, true);
+ gf->model()->writeMSH("multiscalePARTS.msh", 2.2, false, true);
//Remesh new faces (Compound Lines and Compound Surfaces)
//-----------------------------------------------------
@@ -1418,12 +1418,12 @@ void partitionAndRemesh(GFaceCompound *gf)
Msg::Info("*** Starting Mesh of surface %d ...", gf->tag());
- //lloydAlgorithm lloyd(10,0);
+ //lloydAlgorithm
for (int i=0; i < NF; i++){
GFace *gfc = gf->model()->getFaceByTag(numf + NF + i );
meshGFace mgf;
mgf(gfc);
- // gfc->lloyd(20,0);
+ //gfc->lloyd(20,0);
for(unsigned int j = 0; j < gfc->triangles.size(); ++j){
MTriangle *t = gfc->triangles[j];
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index 0ea08b12a2..610dde401f 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -218,6 +218,7 @@ static int getAspectRatio(std::vector<MElement *> &elements,
static void getGenusAndRatio(std::vector<MElement *> &elements, int & genus, int &AR, int &NB)
{
std::vector<std::vector<MEdge> > boundaries;
+ boundaries.clear();
genus = getGenus(elements, boundaries);
NB = boundaries.size();
AR = getAspectRatio(elements, boundaries);
@@ -252,7 +253,7 @@ static void printLevel(std::vector<MElement *> &elements, int recur, int region)
{
char fn[256];
sprintf(fn, "part_%d_%d.msh", recur, region);
- double version = 2.0;
+ double version = 2.2;
std::set<MVertex*> vs;
for (unsigned int i = 0; i < elements.size(); i++){
@@ -293,7 +294,7 @@ multiscalePartition::multiscalePartition(std::vector<MElement *> &elements,
options.num_partitions = nbParts;
options.partitioner = 1; //1 CHACO, 2 METIS
if (options.partitioner == 1){
- options.global_method = 2;// 1 Multilevel-KL, 2 Spectral
+ options.global_method = 1;// 1 Multilevel-KL, 2 Spectral
options.mesh_dims[0] = nbParts;
}
@@ -346,7 +347,7 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts,
int genus, AR, NB;
getGenusAndRatio(regions[i], genus, AR, NB);
- //printLevel (nextLevel->elements, nextLevel->recur,nextLevel->region);
+ printLevel (nextLevel->elements, nextLevel->recur,nextLevel->region);
if (genus < 0) {
Msg::Error("Genus partition is negative G=%d!", genus);
@@ -359,7 +360,7 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts,
nextLevel->recur,nextLevel->region, genus, AR, nbParts);
partition(*nextLevel, nbParts, MULTILEVEL);
}
- else if (genus == 0 && AR > 3 || genus == 0 && NB > 1){
+ else if (genus == 0 && AR > 5 ){// || genus == 0 && NB > 1){
int nbParts = 2;
Msg::Info("Mesh partition: level (%d-%d) is ZERO-GENUS (AR=%d NB=%d) ---> LAPLACIAN partition %d parts",
nextLevel->recur,nextLevel->region, AR, NB, nbParts);
diff --git a/Solver/TESTCASES/Aorta.lua b/Solver/TESTCASES/Aorta.lua
new file mode 100644
index 0000000000..4f475146d6
--- /dev/null
+++ b/Solver/TESTCASES/Aorta.lua
@@ -0,0 +1,67 @@
+--[[
+ Function for initial conditions
+--]]
+
+
+function initial_condition( xyz , f )
+ for i=0,xyz:size1()-1 do
+ f:set (i, 0, 3.14)
+ f:set (i, 1, 0)
+ end
+end
+
+PI = 3.14159
+T= 0.05
+Amax = 3.15
+
+function inlet( FCT )
+ FCT:set(0,0,Amax*math.sin(2*PI*t/T)*math.step(T/2-t))
+ FCT:set(0,1, 0)
+end
+
+xyz = functionCoordinates.get()
+
+--[[
+ Example of a lua program driving the DG code
+--]]
+
+model = GModel()
+model:load ('aorta.msh')
+order=3
+dimension=1
+
+-- boundary condition
+inlet_bc = functionLua(2, 'inlet')
+law = dgConservationLawShallowWater1d()
+law:addBoundaryCondition('Inlet',law:newOutsideValueBoundary(inlet_bc))
+law:addBoundaryCondition('Outlet',law:newBoundaryWall())
+law:setBathymetry(functionConstant({0}))
+law:setLinearDissipation(functionConstant({0}))
+
+groups = dgGroupCollection(model, dimension, order)
+groups:buildGroupsOfInterfaces()
+
+rk=dgRungeKutta()
+limiter = dgSlopeLimiter(law)
+
+-- build solution vector
+FS = functionLua(2, 'initial_condition', {xyz})
+solution = dgDofContainer(groups, law:getNbFields())
+solution:L2Projection(FS)
+
+print'*** print initial sol ***'
+solution:exportMsh('output/init')
+
+print'*** solve ***'
+--dt = 0.00001;
+CFL = 0.2;
+for i=1,9000 do
+ dt = CFL * rk:computeInvSpectralRadius(law,solution);
+ norm = rk:iterate44(law,dt,solution)
+ if (i % 100 == 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
diff --git a/Solver/TESTCASES/CylinderEddies.lua b/Solver/TESTCASES/CylinderEddies.lua
index 33c66654ef..57347cdd5a 100644
--- a/Solver/TESTCASES/CylinderEddies.lua
+++ b/Solver/TESTCASES/CylinderEddies.lua
@@ -21,8 +21,8 @@ end
--]]
order = 5
-
-FS = functionLua(4, 'free_stream', {'XYZ'}):getName()
+xyz = functionCoordinates.get()
+FS = functionLua(4, 'free_stream', {xyz})
-- diffusivity
mu=fullMatrix(1,1);
diff --git a/Solver/TESTCASES/DamBreak.lua b/Solver/TESTCASES/DamBreak.lua
index a192e7ba45..2c5b1da992 100644
--- a/Solver/TESTCASES/DamBreak.lua
+++ b/Solver/TESTCASES/DamBreak.lua
@@ -30,39 +30,38 @@ dimension=1
law = dgConservationLawShallowWater1d()
law:addBoundaryCondition('Left',law:newBoundaryWall())
law:addBoundaryCondition('Right',law:newBoundaryWall())
+law:setBathymetry(functionConstant({0}))
+law:setLinearDissipation(functionConstant({0}))
groups = dgGroupCollection(model, dimension, order)
groups:buildGroupsOfInterfaces()
-claw = dgConservationLawShallowWater1d()
-claw:addBoundaryCondition('Wall',law:newBoundaryWall())
-
rk=dgRungeKutta()
limiter = dgSlopeLimiter(law)
rk:setLimiter(limiter)
-- build solution vector
xyz = functionCoordinates.get()
-FS = functionLua(1, 'initial_condition', {xyz})
+FS = functionLua(2, 'initial_condition', {xyz})
solution = dgDofContainer(groups, law:getNbFields())
solution:L2Projection(FS)
print'*** print initial sol ***'
-solution:exportMsh('output/init')
+--solution:exportMsh('output/init')
limiter:apply(solution)
solution:exportMsh('output/init_limit')
print'*** solve ***'
--dt = 0.00001;
-CFL = 1;
-for i=1,1000 do
+CFL = 0.2;
+for i=1,15000 do
dt = CFL * rk:computeInvSpectralRadius(law,solution);
--- norm = rk:iterate44(law,dt,solution)
- norm = rk:iterateEuler(law,dt,solution)
- if (i % 20 == 0) then
+ norm = rk:iterate44(law,dt,solution)
+-- norm = rk:iterateEuler(law,dt,solution)
+ if (i % 100 == 0) then
print('|ITER|',i,'|NORM|',norm,'|DT|',dt,'|CPU|',os.clock() - x)
end
- if (i % 20 == 0) then
+ if (i % 100 == 0) then
solution:exportMsh(string.format('output/solution-%06d', i))
end
end
diff --git a/Solver/TESTCASES/DamBreak2d.lua b/Solver/TESTCASES/DamBreak2d.lua
new file mode 100644
index 0000000000..876827531e
--- /dev/null
+++ b/Solver/TESTCASES/DamBreak2d.lua
@@ -0,0 +1,72 @@
+--[[
+ 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.0) then
+ f:set (i, 0, 40)
+ f:set (i, 1, 0)
+ f:set (i, 2, 0)
+ else
+ f:set (i, 0, 5)
+ f:set (i, 1, 0)
+ f:set (i, 2, 0)
+ end
+ end
+end
+
+--[[
+ Example of a lua program driving the DG code
+--]]
+
+model = GModel()
+model:load ('rect2.msh')
+order=1
+dimension=2
+
+-- boundary condition
+law = dgConservationLawShallowWater2d()
+law:addBoundaryCondition('Wall',law:newBoundaryWall())
+law:addBoundaryCondition('Symmetry',law:newSymmetryBoundary())
+law:setCoriolisFactor(functionConstant({0}))
+law:setQuadraticDissipation(functionConstant({0}))
+law:setLinearDissipation(functionConstant({0}))
+law:setSource(functionConstant({0}))
+law:setBathymetry(functionConstant({0}))
+
+groups = dgGroupCollection(model, dimension, order)
+groups:buildGroupsOfInterfaces()
+
+rk=dgRungeKutta()
+limiter = dgSlopeLimiter(law)
+rk:setLimiter(limiter)
+
+-- build solution vector
+xyz = functionCoordinates.get()
+FS = functionLua(1, 'initial_condition', {xyz})
+solution = dgDofContainer(groups, law:getNbFields())
+solution:L2Projection(FS)
+
+print'*** print initial sol ***'
+--solution:exportMsh('output/init')
+limiter:apply(solution)
+solution:exportMsh('output/init_limit')
+
+print'*** solve ***'
+--dt = 0.00001;
+CFL = 0.02;
+for i=1,9000 do
+ dt = CFL * rk:computeInvSpectralRadius(law,solution);
+-- norm = rk:iterate44(law,dt,solution)
+ norm = rk:iterateEuler(law,dt,solution)
+ if (i % 100 == 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
diff --git a/Solver/TESTCASES/ForwardFacingStepNEW.lua b/Solver/TESTCASES/ForwardFacingStepNEW.lua
index 1eb870f462..a87e559bfb 100644
--- a/Solver/TESTCASES/ForwardFacingStepNEW.lua
+++ b/Solver/TESTCASES/ForwardFacingStepNEW.lua
@@ -35,7 +35,7 @@ model:load ('step.msh')
order = 1
dimension = 2
-FS = functionLua(4, 'free_stream', {'XYZ'}):getName()
+FS = functionLua(4, 'free_stream', {XYZ})
-- boundary condition
law=dgPerfectGasLaw2d()
diff --git a/Solver/TESTCASES/aorta.geo b/Solver/TESTCASES/aorta.geo
new file mode 100644
index 0000000000..70ea0e186a
--- /dev/null
+++ b/Solver/TESTCASES/aorta.geo
@@ -0,0 +1,8 @@
+
+Point(1) = {0, 0, 0, 0.5};
+Point(2) = {100, 0, 0, 0.5};
+Line(1) = {1, 2};
+
+Physical Point("Inlet") = {1};
+Physical Point("Outlet") = {2};
+Physical Line("Line") = {1};
diff --git a/Solver/TESTCASES/edge.msh b/Solver/TESTCASES/edge.msh
index 2f056713bd..cc22a33bec 100644
--- a/Solver/TESTCASES/edge.msh
+++ b/Solver/TESTCASES/edge.msh
@@ -8,30 +8,210 @@ $PhysicalNames
1 3 "Line"
$EndPhysicalNames
$Nodes
-20
+200
1 -1 0 0
2 1 0 0
-3 -0.8947368421055554 0 0
-4 -0.7894736842111106 0 0
-5 -0.684210526316666 0 0
-6 -0.5789473684222213 0 0
-7 -0.4736842105277401 0 0
-8 -0.3684210526331494 0 0
-9 -0.2631578947385587 0 0
-10 -0.1578947368439678 0 0
-11 -0.0526315789493772 0 0
-12 0.05263157894535975 0 0
-13 0.1578947368402426 0 0
-14 0.2631578947351256 0 0
-15 0.3684210526300085 0 0
-16 0.4736842105248913 0 0
-17 0.5789473684198838 0 0
-18 0.6842105263149127 0 0
-19 0.7894736842099421 0 0
-20 0.8947368421049708 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.8693467336687043 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.8090452261311833 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.7386934673374086 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.6783919597998876 0 0
+35 -0.6683417085436341 0 0
+36 -0.6582914572873806 0 0
+37 -0.6482412060311271 0 0
+38 -0.6381909547748735 0 0
+39 -0.6281407035186201 0 0
+40 -0.6180904522623666 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.5577889447248455 0 0
+47 -0.5477386934685919 0 0
+48 -0.5376884422123385 0 0
+49 -0.527638190956085 0 0
+50 -0.5175879396998314 0 0
+51 -0.5075376884435779 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.3366834170870416 0 0
+69 -0.3266331658307741 0 0
+70 -0.3165829145745066 0 0
+71 -0.3065326633182393 0 0
+72 -0.2964824120619718 0 0
+73 -0.2864321608057042 0 0
+74 -0.2763819095494368 0 0
+75 -0.2663316582931694 0 0
+76 -0.2562814070369018 0 0
+77 -0.2462311557806344 0 0
+78 -0.236180904524367 0 0
+79 -0.2261306532680996 0 0
+80 -0.216080402011832 0 0
+81 -0.2060301507555646 0 0
+82 -0.1959798994992972 0 0
+83 -0.1859296482430296 0 0
+84 -0.1758793969867622 0 0
+85 -0.1658291457304948 0 0
+86 -0.1557788944742273 0 0
+87 -0.1457286432179599 0 0
+88 -0.1356783919616924 0 0
+89 -0.1256281407054249 0 0
+90 -0.1155778894491575 0 0
+91 -0.10552763819289 0 0
+92 -0.09547738693662255 0 0
+93 -0.08542713568035498 0 0
+94 -0.07537688442408763 0 0
+95 -0.06532663316782017 0 0
+96 -0.05527638191155271 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.01507537688236127 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.07537688442013368 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.175879396983087 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.2763819095460405 0 0
+130 0.2864321608023359 0 0
+131 0.2964824120586314 0 0
+132 0.3065326633149268 0 0
+133 0.3165829145712222 0 0
+134 0.3266331658275172 0 0
+135 0.3366834170838127 0 0
+136 0.3467336683401081 0 0
+137 0.3567839195964035 0 0
+138 0.366834170852699 0 0
+139 0.3768844221089944 0 0
+140 0.3869346733652894 0 0
+141 0.3969849246215849 0 0
+142 0.4070351758778803 0 0
+143 0.4170854271341757 0 0
+144 0.427135678390471 0 0
+145 0.4371859296467664 0 0
+146 0.4472361809030618 0 0
+147 0.457286432159357 0 0
+148 0.4673366834156523 0 0
+149 0.4773869346719477 0 0
+150 0.4874371859282431 0 0
+151 0.4974874371845386 0 0
+152 0.5075376884408445 0 0
+153 0.5175879396971534 0 0
+154 0.5276381909534629 0 0
+155 0.5376884422097721 0 0
+156 0.5477386934660815 0 0
+157 0.5577889447223909 0 0
+158 0.5678391959787001 0 0
+159 0.5778894472350093 0 0
+160 0.5879396984913186 0 0
+161 0.597989949747628 0 0
+162 0.6080402010039372 0 0
+163 0.6180904522602466 0 0
+164 0.628140703516556 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.6984924623107212 0 0
+172 0.7085427135670301 0 0
+173 0.7185929648233396 0 0
+174 0.7286432160796488 0 0
+175 0.7386934673359582 0 0
+176 0.7487437185922676 0 0
+177 0.7587939698485768 0 0
+178 0.768844221104886 0 0
+179 0.7788944723611952 0 0
+180 0.7889447236175047 0 0
+181 0.7989949748738141 0 0
+182 0.8090452261301233 0 0
+183 0.8190954773864327 0 0
+184 0.8291457286427419 0 0
+185 0.8391959798990509 0 0
+186 0.8492462311553604 0 0
+187 0.8592964824116698 0 0
+188 0.869346733667979 0 0
+189 0.8793969849242884 0 0
+190 0.8894472361805978 0 0
+191 0.8994974874369068 0 0
+192 0.9095477386932163 0 0
+193 0.9195979899495255 0 0
+194 0.9296482412058349 0 0
+195 0.9396984924621443 0 0
+196 0.9497487437184535 0 0
+197 0.959798994974763 0 0
+198 0.9698492462310719 0 0
+199 0.9798994974873814 0 0
+200 0.9899497487436908 0 0
$EndNodes
$Elements
-21
+201
1 15 2 1 1 1
2 15 2 2 2 2
3 1 2 3 1 1 3
@@ -52,5 +232,185 @@ $Elements
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 2
+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/dgConservationLawShallowWater1d.cpp b/Solver/dgConservationLawShallowWater1d.cpp
index f56a1dd1b6..f24865f417 100644
--- a/Solver/dgConservationLawShallowWater1d.cpp
+++ b/Solver/dgConservationLawShallowWater1d.cpp
@@ -1,19 +1,87 @@
#include "dgConservationLawShallowWater1d.h"
#include "function.h"
static double g = 9.81;
-static double h = 0.;
-static double linearDissipation = 0.0;
+
+class dgConservationLawShallowWater1d : public dgConservationLaw {
+ class advection;
+ class source;
+ class riemann;
+ class boundaryWall;
+ class clipToPhysics;
+ class maxConvectiveSpeed;
+ const function *_bathymetry, *_linearDissipation;
+ public:
+ dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const;
+ dataCacheDouble *newMaxConvectiveSpeed( dataCacheMap &cacheMap) const;
+ dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const;
+ dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const;
+ dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const;
+ dataCacheDouble *newClipToPhysics (dataCacheMap &cacheMap) const;
+ inline void setBathymetry(const function *bathymetry) {_bathymetry = bathymetry;}
+ inline void setLinearDissipation(const function *linearDissipation){_linearDissipation = linearDissipation;}
+ dgConservationLawShallowWater1d()
+ {
+ _nbf = 2; // eta u
+ }
+ dgBoundaryCondition *newBoundaryWall();
+};
+
+class dgConservationLawShallowWater1d::clipToPhysics : public dataCacheDouble {
+ dataCacheDouble &sol, &_bathymetry;
+ double _hMin;
+public:
+ clipToPhysics(dataCacheMap &cacheMap, const function *bathymetry, double hMin):
+ dataCacheDouble(cacheMap,1,3),
+ sol(cacheMap.getSolution(this)),
+ _bathymetry(cacheMap.get(bathymetry,this))
+ {
+ _hMin=hMin;
+ };
+ void _eval () {
+ int nQP = _value.size1();
+ for (size_t k = 0 ; k < nQP; k++ ){
+ double h = _bathymetry(k,0);
+ _value(k,0) = sol(k,0);
+ _value(k,1) = sol(k,1);
+ double H = sol(k,0)+h;
+ if (H < _hMin){
+ _value(k,0) = _hMin;
+ }
+ }
+ }
+};
+
+class dgConservationLawShallowWater1d::maxConvectiveSpeed : public dataCacheDouble {
+ dataCacheDouble &sol, &_bathymetry;
+ public:
+ maxConvectiveSpeed(dataCacheMap &cacheMap, const function *bathymetry):
+ dataCacheDouble(cacheMap,1,1),
+ sol(cacheMap.getSolution(this)),
+ _bathymetry(cacheMap.get(bathymetry,this))
+ {};
+ void _eval () {
+ int nQP = sol().size1();
+ for(int i=0; i< nQP; i++) {
+ double h = _bathymetry(i,0);
+ double eta = sol(i,0);
+ const double c = sqrt(g*eta);
+ _value(i,0) = sol(i,1) + c;
+ }
+ }
+};
class dgConservationLawShallowWater1d::advection: public dataCacheDouble {
- dataCacheDouble /
+ dataCacheDouble &sol, &_bathymetry;
public:
- advection(dataCacheMap &cacheMap):
+ advection(dataCacheMap &cacheMap, const function *bathymetry):
dataCacheDouble(cacheMap,1,6),
- sol(cacheMap.getSolution(this))
+ sol(cacheMap.getSolution(this)),
+ _bathymetry(cacheMap.get(bathymetry,this))
{};
void _eval () {
int nQP = _value.size1();
for(int i=0; i< nQP; i++) {
+ double h = _bathymetry(i,0);
double eta = sol(i,0);
double u = sol(i,1);
// flux_x
@@ -29,54 +97,41 @@ class dgConservationLawShallowWater1d::advection: public dataCacheDouble {
}
};
-class dgConservationLawShallowWater1d::maxConvectiveSpeed : public dataCacheDouble {
- dataCacheDouble /
- public:
- maxConvectiveSpeed(dataCacheMap &cacheMap):
- dataCacheDouble(cacheMap,1,1),
- sol(cacheMap.getSolution(this))
- {
- };
- void _eval () {
- int nQP = sol().size1();
- for(int i=0; i< nQP; i++) {
- const double c = sqrt(g*sol(i,0));
- _value(i,0) = sol(i,1) + c;
- }
- }
-};
-
class dgConservationLawShallowWater1d::source: public dataCacheDouble {
- dataCacheDouble &xyz, &solution,&solutionGradient;
+ dataCacheDouble &solution;
+ dataCacheDouble &_linearDissipation, &_bathymetry;
public :
- source(dataCacheMap &cacheMap) :
+ source(dataCacheMap &cacheMap, const function *linearDissipation, const function *bathymetry) :
dataCacheDouble(cacheMap,1,2),
- xyz(cacheMap.getParametricCoordinates(this)),
+ _linearDissipation(cacheMap.get(linearDissipation,this)),
solution(cacheMap.getSolution(this)),
- solutionGradient(cacheMap.getSolutionGradient(this))
+ _bathymetry(cacheMap.get(bathymetry,this))
{}
void _eval () {
int nQP =_value.size1();
for (int i=0; i<nQP; i++) {
+ double h = _bathymetry(i,0);
double eta = solution(i,0);
double u = solution(i,1);
_value (i,0) = 0;
- _value (i,1) = - linearDissipation*u ;
+ _value (i,1) = - _linearDissipation(i,0)*u/(h+eta) ;
}
}
};
class dgConservationLawShallowWater1d::riemann:public dataCacheDouble {
- dataCacheDouble &normals, &solL, &solR;
+ dataCacheDouble &normals, &solL, &solR, &_bathymetry;;
public:
- riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
+ riemann(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight, const function *bathymetry):
dataCacheDouble(cacheMapLeft,1,4),
normals(cacheMapLeft.getNormals(this)),
solL(cacheMapLeft.getSolution(this)),
- solR(cacheMapRight.getSolution(this))
+ solR(cacheMapRight.getSolution(this)),
+ _bathymetry(cacheMapLeft.get(bathymetry, this))
{};
void _eval () {
int nQP = solL().size1();
for(int i=0; i< nQP; i++) {
+ double h = _bathymetry(i,0);
double nx = normals(0,i);
double HL = solL(i,0) + h, HR = solR(i,0) + h;
double uL = nx*solL(i,1) , uR = nx*solR(i,1) ;
@@ -129,44 +184,20 @@ class dgConservationLawShallowWater1d::boundaryWall : public dgBoundaryCondition
};
-
-class dgConservationLawShallowWater1d::clipToPhysics : public dataCacheDouble {
- dataCacheDouble /
- double _hMin;
-public:
- clipToPhysics(dataCacheMap &cacheMap, double hMin):
- dataCacheDouble(cacheMap,1,2),
- sol(cacheMap.getSolution(this))
- {
- _hMin=hMin;
- };
- void _eval () {
- const int nDofs = sol().size1();
- for (size_t k = 0 ; k < nDofs; k++ ){
- _value(k,0) = sol(k,0);
- _value(k,1) = sol(k,1);
- double H = sol(k,0)+h;
- if (H < _hMin){
- _value(k,0) = _hMin;
- }
- }
- }
-};
-
dataCacheDouble *dgConservationLawShallowWater1d::newMaxConvectiveSpeed( dataCacheMap &cacheMap) const {
- return new maxConvectiveSpeed(cacheMap);
+ return new maxConvectiveSpeed(cacheMap, _bathymetry);
}
dataCacheDouble *dgConservationLawShallowWater1d::newConvectiveFlux( dataCacheMap &cacheMap) const {
- return new advection(cacheMap);
+ return new advection(cacheMap, _bathymetry);
}
dataCacheDouble *dgConservationLawShallowWater1d::newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {
- return new riemann(cacheMapLeft, cacheMapRight);
+ return new riemann(cacheMapLeft, cacheMapRight, _bathymetry);
}
dataCacheDouble *dgConservationLawShallowWater1d::newDiffusiveFlux( dataCacheMap &cacheMap) const {
return 0;
}
dataCacheDouble *dgConservationLawShallowWater1d::newSourceTerm (dataCacheMap &cacheMap) const {
- return new source(cacheMap);
+ return new source(cacheMap, _linearDissipation, _bathymetry);
}
dgBoundaryCondition *dgConservationLawShallowWater1d::newBoundaryWall(){
@@ -174,7 +205,7 @@ dgBoundaryCondition *dgConservationLawShallowWater1d::newBoundaryWall(){
}
dataCacheDouble *dgConservationLawShallowWater1d::newClipToPhysics( dataCacheMap &cacheMap) const {
- return new clipToPhysics(cacheMap,1e-5);
+ return new clipToPhysics(cacheMap, _bathymetry, 1e-5);
}
#include "Bindings.h"
@@ -187,4 +218,10 @@ void dgConservationLawShallowWater1dRegisterBindings (binding *b){
cm = cb->setConstructor<dgConservationLawShallowWater1d>();
cm->setDescription("A new shallow water 1d conservation law.");
cb->setParentClass<dgConservationLaw>();
+ cm = cb->addMethod("setLinearDissipation", &dgConservationLawShallowWater1d::setLinearDissipation);
+ cm->setDescription("set the function to evaluate the linear dissipation Kr u/A");
+ cm->setArgNames("Kr",NULL);
+ cm = cb->addMethod("setBathymetry", &dgConservationLawShallowWater1d::setBathymetry);
+ cm->setDescription("set the function to evaluate the bathymetry h (H = h+eta)");
+ cm->setArgNames("h",NULL);
}
diff --git a/Solver/dgConservationLawShallowWater1d.h b/Solver/dgConservationLawShallowWater1d.h
index 69690da7fb..a2642a642e 100644
--- a/Solver/dgConservationLawShallowWater1d.h
+++ b/Solver/dgConservationLawShallowWater1d.h
@@ -2,28 +2,9 @@
#define _DG_CONSERVATION_LAW_SHALLOW_WATER_1D_
#include "dgConservationLaw.h"
-class dataCacheMap;
+
class binding;
-class dgConservationLawShallowWater1d : public dgConservationLaw {
- class advection;
- class source;
- class riemann;
- class boundaryWall;
- class clipToPhysics;
- class maxConvectiveSpeed;
- public:
- dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const;
- dataCacheDouble *newMaxConvectiveSpeed( dataCacheMap &cacheMap) const;
- dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const;
- dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const;
- dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const;
- dataCacheDouble *newClipToPhysics (dataCacheMap &cacheMap) const;
- dgConservationLawShallowWater1d()
- {
- _nbf = 2; // eta u
- }
- dgBoundaryCondition *newBoundaryWall();
-};
void dgConservationLawShallowWater1dRegisterBindings(binding *b);
+
#endif
diff --git a/Solver/dgConservationLawShallowWater2d.cpp b/Solver/dgConservationLawShallowWater2d.cpp
index b09a5c69aa..6394f67313 100644
--- a/Solver/dgConservationLawShallowWater2d.cpp
+++ b/Solver/dgConservationLawShallowWater2d.cpp
@@ -6,6 +6,7 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
class advection;
class source;
class riemann;
+ class clipToPhysics;
class boundaryWall;
class maxConvectiveSpeed;
const function *_linearDissipation, *_quadraticDissipation, *_source, *_coriolisFactor, *_bathymetry;
@@ -15,6 +16,7 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const;
dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const;
dataCacheDouble *newMaxConvectiveSpeed (dataCacheMap &cacheMap) const;
+ dataCacheDouble *newClipToPhysics (dataCacheMap &cacheMap) const;
inline void setCoriolisFactor(const function *coriolisFactor){_coriolisFactor = coriolisFactor;}
inline void setLinearDissipation(const function *linearDissipation){_linearDissipation = linearDissipation;}
inline void setQuadraticDissipation(const function *quadraticDissipation){_quadraticDissipation = quadraticDissipation;}
@@ -27,6 +29,32 @@ class dgConservationLawShallowWater2d : public dgConservationLaw {
dgBoundaryCondition *newBoundaryWall();
};
+class dgConservationLawShallowWater2d::clipToPhysics : public dataCacheDouble {
+ dataCacheDouble &sol, &_bathymetry;
+ double _hMin;
+public:
+ clipToPhysics(dataCacheMap &cacheMap, const function *bathymetry, double hMin):
+ dataCacheDouble(cacheMap,1,3),
+ sol(cacheMap.getSolution(this)),
+ _bathymetry(cacheMap.get(bathymetry,this))
+ {
+ _hMin=hMin;
+ };
+ void _eval () {
+ const int nQP = _value.size1();
+ for (size_t k = 0 ; k < nQP; k++ ){
+ double h = _bathymetry(k,0);
+ _value(k,0) = sol(k,0);
+ _value(k,1) = sol(k,1);
+ _value(k,2) = sol(k,2);
+ double H = sol(k,0)+h;
+ if (H < _hMin){
+ _value(k,0) = _hMin;
+ }
+ }
+ }
+};
+
class dgConservationLawShallowWater2d::maxConvectiveSpeed: public dataCacheDouble {
dataCacheDouble &sol, &_bathymetry;
public:
@@ -227,6 +255,10 @@ dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryWall(){
return new boundaryWall(this);
}
+dataCacheDouble *dgConservationLawShallowWater2d::newClipToPhysics( dataCacheMap &cacheMap) const {
+ return new clipToPhysics(cacheMap, _bathymetry, 1e-5);
+}
+
#include "Bindings.h"
void dgConservationLawShallowWater2dRegisterBindings (binding *b){
classBinding *cb = b->addClass<dgConservationLawShallowWater2d>("dgConservationLawShallowWater2d");
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 1b7d547b28..adbfdf1b2e 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -249,8 +249,7 @@ static void recur_compute_centers_ (double R, double a1, double a2,
SPoint2 PL (R*cos(a1),R*sin(a1));
SPoint2 PR (R*cos(a2),R*sin(a2));
- centers.push_back(std::make_pair(PL,zero));
- centers.push_back(std::make_pair(PR,zero));
+
std::vector<SPoint2> centersChild;
centersChild.clear();
@@ -271,8 +270,9 @@ static void recur_compute_centers_ (double R, double a1, double a2,
//add the center of real holes ...
std::vector<std::vector<MEdge> > boundaries;
connected_bounds(root->elements, boundaries);
- int toadd = 0;
- if (root->children.size()==0 || boundaries.size()-1 != root->children.size() ){
+ int added = 0;
+ int toadd = boundaries.size()-1 - root->children.size();
+ if (root->children.size()==0 || (root->children.size()> 0 && toadd > 0) ){
for (unsigned int i = 0; i < boundaries.size(); i++){
std::vector<MEdge> me = boundaries[i];
SPoint2 c(0.0,0.0);
@@ -295,22 +295,24 @@ static void recur_compute_centers_ (double R, double a1, double a2,
SPoint2 p = *it2;
double dist = sqrt ((c.x() - p.x())*(c.x() - p.x())+
(c.y() - p.y())*(c.y() - p.y()));
- if (dist < 0.6*rad) newCenter = false;//0.6
+ if (dist < 0.5*rad) newCenter = false;//0.6
}
- if (std::abs(rad/root->radius) < 0.6 && std::abs(rad) < 0.95 && newCenter){//0.6
- toadd++;
+ if (std::abs(rad/root->radius) < 0.65 && std::abs(rad) < 0.95 && newCenter){//0.6
+ added++;
centers.push_back(std::make_pair(c,zero));
}
}
}
- if (toadd != boundaries.size()-1-root->children.size() ) {
- printf("!!!!!!!! ARG added =%d != %d\n", toadd, boundaries.size()-1- root->children.size());
+ if (added != toadd && root->children.size()> 0) {
+ printf("!!!!!!!! ARG added =%d != %d (bounds=%d, child=%d)\n", added, boundaries.size()-1- root->children.size(), boundaries.size(), root->children.size());
//exit(1);
}
//sort centers
std::sort(centers.begin(),centers.end(), sort_pred(PL,PR));
+ centers.insert(centers.begin(), std::make_pair(PL,zero));
+ centers.push_back(std::make_pair(PR,zero));
for (int i=1;i<centers.size()-1;i++){
multiscaleLaplaceLevel* m1 = centers[i-1].second;
@@ -820,7 +822,7 @@ multiscaleLaplace::multiscaleLaplace (std::vector<MElement *> &elements,
//Compute centers for the cut
int nbElems = 0;
recur_compute_centers_ (1.0, M_PI, 0.0, root, nbElems);
- printf("CENTERS: elements =%d, recur nbElems = %d \n", elements.size(), nbElems);
+ //printf("CENTERS: elements =%d, recur nbElems = %d \n", elements.size(), nbElems);
//Partition the mesh in left and right
cut (elements);
@@ -943,7 +945,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) //1.e-5
+ if (local_size < 5.e-4*global_size) //1.e-5
tooSmall.push_back(e);
else goodSize.push_back(e);
}
@@ -992,7 +994,7 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
// double area = fabs(triangle_area(q0, q1, q2));
// totArea += area;
double local_size = localSize(e,solution);
- if (local_size < 5.e-7 * global_size) //1.e-7
+ if (local_size < 5.e-6 * global_size) //1.e-7
really_small_elements = true;
}
//center *= (1./regions_[i].size());
@@ -1174,6 +1176,10 @@ void multiscaleLaplace::cut (std::vector<MElement *> &elements)
recur_cut_ (1.0, M_PI, 0.0, root,left,right);
connected_left_right(left, right);
+ printLevel ("Rootcut-left.msh",left,0,2.2);
+ printLevel ("Rootcut-right.msh",right,0,2.2);
+ printLevel ("Rootcut-all.msh",elements, 0,2.2);
+
if ( elements.size() != left.size()+right.size()) {
Msg::Error("Cutting laplace wrong nb elements (%d) != left + right (%d)", elements.size(), left.size()+right.size());
exit(1);
@@ -1183,9 +1189,6 @@ void multiscaleLaplace::cut (std::vector<MElement *> &elements)
elements.insert(elements.end(),left.begin(),left.end());
elements.insert(elements.end(),right.begin(),right.end());
- printLevel ("Rootcut-left.msh",left,0,2.2);
- printLevel ("Rootcut-right.msh",right,0,2.2);
- printLevel ("Rootcut-all.msh",elements, 0,2.2);
//exit(1);
}
diff --git a/benchmarks/3d/Torus.geo b/benchmarks/3d/Torus.geo
index 46241c987f..56489b2a6a 100644
--- a/benchmarks/3d/Torus.geo
+++ b/benchmarks/3d/Torus.geo
@@ -1,4 +1,4 @@
-lc = .9;
+lc = 15.5;
Point(1) = {2.0,0.0,0.0,lc};
Point(2) = {2.0,1,0.0,lc};
Point(3) = {1,0,0.0,lc};
@@ -12,5 +12,5 @@ Circle(4) = {5,1,4};
Line Loop(5) = {4,1,2,3};
Plane Surface(6) = {5};
-Extrude Surface{6, {0.0,1,0}, {0,0.0,0.0}, 1*3.14159/2};
+Extrude Surface{6, {0.0,1,0}, {0,0.0,0.0}, 1*3.14159};
//Recombine Surface {6, 27, 23, 15, 19, 28};
--
GitLab