diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index f3d55e39149840a5c395bfb2264340a049ecc2f9..9b43f6340ca65a7404507613a24821c9d877e889 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -752,23 +752,24 @@ static void generate3dFaceClosure(polynomialBasis::clCont &closure, int order)
}
}
-static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order)
+static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order, int nNod = 3)
{
closure.clear();
- closure.resize(6);
- for (int j = 0; j < 3 ; j++){
+ closure.resize(2*nNod);
+ for (int j = 0; j < nNod ; j++){
closure[j].push_back(j);
- closure[j].push_back((j+1)%3);
- closure[3+j].push_back((j+1)%3);
- closure[3+j].push_back(j);
+ closure[j].push_back((j+1)%nNod);
+ closure[nNod+j].push_back((j+1)%nNod);
+ closure[nNod+j].push_back(j);
for (int i=0; i < order-1; i++){
- closure[j].push_back( 3 + (order-1)*j + i );
- closure[3+j].push_back(3 + (order-1)*(j+1) -i -1);
+ closure[j].push_back( nNod + (order-1)*j + i );
+ closure[nNod+j].push_back(nNod + (order-1)*(j+1) -i -1);
}
}
}
+
static void generate1dVertexClosure(polynomialBasis::clCont &closure)
{
closure.clear();
@@ -893,38 +894,47 @@ const polynomialBasis &polynomialBases::find(int tag)
case MSH_QUA_4 :
F.monomials = generatePascalQuad(1);
F.points = gmshGeneratePointsQuad(1,false);
+ generate2dEdgeClosure(F.edgeClosure, 1, 4);
break;
case MSH_QUA_9 :
F.monomials = generatePascalQuad(2);
F.points = gmshGeneratePointsQuad(2,false);
+ generate2dEdgeClosure(F.edgeClosure, 2, 4);
break;
case MSH_QUA_16 :
F.monomials = generatePascalQuad(3);
F.points = gmshGeneratePointsQuad(3,false);
+ generate2dEdgeClosure(F.edgeClosure, 3, 4);
break;
case MSH_QUA_25 :
F.monomials = generatePascalQuad(4);
F.points = gmshGeneratePointsQuad(4,false);
+ generate2dEdgeClosure(F.edgeClosure, 4, 4);
break;
case MSH_QUA_36 :
F.monomials = generatePascalQuad(5);
F.points = gmshGeneratePointsQuad(5,false);
+ generate2dEdgeClosure(F.edgeClosure, 5, 4);
break;
case MSH_QUA_8 :
F.monomials = generatePascalQuadSerendip(2);
F.points = gmshGeneratePointsQuad(2,true);
+ generate2dEdgeClosure(F.edgeClosure, 2, 4);
break;
case MSH_QUA_12 :
F.monomials = generatePascalQuadSerendip(3);
F.points = gmshGeneratePointsQuad(3,true);
+ generate2dEdgeClosure(F.edgeClosure, 3, 4);
break;
case MSH_QUA_16I :
F.monomials = generatePascalQuadSerendip(4);
F.points = gmshGeneratePointsQuad(4,true);
+ generate2dEdgeClosure(F.edgeClosure, 4, 4);
break;
case MSH_QUA_20 :
F.monomials = generatePascalQuadSerendip(5);
F.points = gmshGeneratePointsQuad(5,true);
+ generate2dEdgeClosure(F.edgeClosure, 5, 4);
break;
default :
Msg::Error("Unknown function space %d: reverting to TET_4", tag);
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index a08c1052feafca865988c411917b513ed016beb4..215bb6cd57eb32529ca963031c0857b5be3b57f8 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -32,7 +32,7 @@ struct polynomialBasis
return faceClosure[id];
}
inline int getEdgeClosureId(int iEdge, int iSign) const {
- return iSign == 1 ? iEdge : 3 + iEdge;
+ return iSign == 1 ? iEdge : edgeClosure.size()/2 + iEdge;
}
inline const std::vector<int> &getEdgeClosure(int id) const
{
diff --git a/Solver/TESTCASES/WavePulse.lua b/Solver/TESTCASES/WavePulse.lua
index 5c679629b912b28c95a275d2e820c6b2dbc5980f..a9251261f9846e87d706631612ae3d440b19f9d0 100644
--- a/Solver/TESTCASES/WavePulse.lua
+++ b/Solver/TESTCASES/WavePulse.lua
@@ -1,3 +1,5 @@
+
+
--[[
Function for initial conditions
--]]
@@ -6,7 +8,11 @@ function initial_condition( xyz , f )
X = xyz:get(i,0) - .5
Y = xyz:get(i,1) - .5
Z = xyz:get(i,2)
- VALUE = math.exp(-40*(X*X+Y*Y+Z*Z));
+-- X2 = xyz:get(i,0) - 1.5
+-- Y2 = xyz:get(i,1) - .5
+-- Z2 = xyz:get(i,2)
+ VALUE = math.exp(-40*(X*X+Y*Y+Z*Z))
+-- + math.exp(-40*(X2*X2+Y2*Y2+Z2*Z2));
f:set(i,0,VALUE)
f:set(i,1,0.0)
f:set(i,2,0.0)
@@ -18,11 +24,12 @@ print'*** Loading the mesh and the model ***'
myModel = GModel ()
--myModel:load('box.geo')
--myModel:load('box.msh')
-myModel:load('square.msh')
+--myModel:load('square_quads.msh')
+myModel:load('square_mixed.msh')
print'*** Create a dg solver ***'
DG = dgSystemOfEquations (myModel)
-DG:setOrder(3)
+DG:setOrder(2)
law=ConservationLawWaveEquation(2)
DG:setConservationLaw(law)
law:addBoundaryCondition('Border',law:newWallBoundary())
@@ -40,12 +47,15 @@ DG:exportSolution('output/solution_000')
print'*** solve ***'
-dt = 0.00125;
N = 1000;
+CFL = 2.1;
+dt = CFL * DG:computeInvSpectralRadius();
+print('DT = ',dt)
for i=1,N do
- norm = DG:multirateRK43(dt)
+-- norm = DG:multirateRK43(dt)
+ norm = DG:RK44(dt)
print('*** ITER ***',i,norm)
- if (i % 10 == 0) then
+ if (i % 100 == 0) then
DG:exportSolution(string.format("output/solution-%04d", i))
end
end
diff --git a/Solver/TESTCASES/square_mixed.geo b/Solver/TESTCASES/square_mixed.geo
new file mode 100644
index 0000000000000000000000000000000000000000..267a1bd2df059225386757395d917c745cc0539c
--- /dev/null
+++ b/Solver/TESTCASES/square_mixed.geo
@@ -0,0 +1,26 @@
+N=25;
+Point(1) = {0.0, 0.0, 0, .03};
+Point(2) = {1, 0.0, 0, .03};
+Point(3) = {1, 1, 0, .03};
+Point(4) = {0, 1, 0, .03};
+Point(5) = {2, 0.0, 0, .03};
+Point(6) = {2, 1, 0, .03};
+Line(1) = {4, 3};
+Line(2) = {3, 2};
+Line(3) = {2, 1};
+Line(4) = {1, 4};
+Line Loop(5) = {2, 3, 4, 1};
+Plane Surface(6) = {5};
+Transfinite Line {1, 2, 4, 3} = N Using Progression 1;
+Transfinite Surface {6};
+Recombine Surface {6};
+//Line(22) = {3, 2};
+Line(7) = {3, 6};
+Line(8) = {6, 5};
+Line(9) = {5, 2};
+Line Loop(10) = {8, 9, -2, 7};
+Plane Surface(11) = {10};
+Transfinite Line {7, 8, 9} = N Using Progression 1;
+Transfinite Surface {11};
+Physical Surface("Inside") = {6,11};
+Physical Line("Border") = {1, 3, 4, 7, 8, 9};
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index a65907abdc4ca3ce80d8101638e0abb4b0f33b09..ff7dcf9c5e046aa284c8e50bd1bca38fe22550aa 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -504,6 +504,22 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
delete boundaryTerm;
}
+/*
+void dgAlgorithm::buildGroupsOfFaces(GModel *model, int dim, int order,
+ std::vector<dgGroupOfElements*> &eGroups,
+ std::vector<dgGroupOfFaces*> &fGroups,
+ std::vector<dgGroupOfFaces*> &bGroups){
+}
+
+void dgAlgorithm::buildMandatoryGroups(dgGroupOfElements &eGroup,
+ std::vector<dgGroupOfElements*> &partitionedGroups){
+}
+
+void dgAlgorithm::partitionGroup(dgGroupOfElements &eGroup,
+ std::vector<dgGroupOfElements*> &partitionedGroups){
+}
+*/
+
void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
std::vector<dgGroupOfElements*> &eGroups,
std::vector<dgGroupOfFaces*> &fGroups,
@@ -514,7 +530,7 @@ void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
std::map<const std::string,std::set<MFace, Less_Face> > boundaryFaces;
std::vector<GEntity*> entities;
model->getEntities(entities);
- std::vector<MElement *> allElements;
+ std::map<int, std::vector<MElement *> >allElements;
for(unsigned int i = 0; i < entities.size(); i++){
GEntity *entity = entities[i];
if(entity->dim() == dim-1){
@@ -539,36 +555,51 @@ void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
}
}else if(entity->dim() == dim){
for (int iel=0; iel<entity->getNumMeshElements(); iel++)
- allElements.push_back(entity->getMeshElement(iel));
+ allElements[entity->getMeshElement(iel)->getType()].push_back(entity->getMeshElement(iel));
}
}
- eGroups.push_back(new dgGroupOfElements(allElements,order));
- switch(dim) {
+ // do a group of elements for every element type that is present in the mesh
+ Msg::Info("%d groups of elements",allElements.size());
+ for (std::map<int, std::vector<MElement *> >::iterator it = allElements.begin(); it !=allElements.end() ; ++it)
+ eGroups.push_back(new dgGroupOfElements(it->second,order));
+
+ for (int i=0;i<eGroups.size();i++){
+ // create a group of faces on the boundary for every group ef elements
+ switch(dim) {
case 1 : {
std::map<const std::string, std::set<MVertex*> >::iterator mapIt;
for(mapIt=boundaryVertices.begin(); mapIt!=boundaryVertices.end(); mapIt++) {
- bGroups.push_back(new dgGroupOfFaces(*eGroups[0],mapIt->first,order,mapIt->second));
+ bGroups.push_back(new dgGroupOfFaces(*eGroups[i],mapIt->first,order,mapIt->second));
}
break;
}
case 2 : {
std::map<const std::string, std::set<MEdge, Less_Edge> >::iterator mapIt;
for(mapIt=boundaryEdges.begin(); mapIt!=boundaryEdges.end(); mapIt++) {
- bGroups.push_back(new dgGroupOfFaces(*eGroups[0],mapIt->first,order,mapIt->second));
+ bGroups.push_back(new dgGroupOfFaces(*eGroups[i],mapIt->first,order,mapIt->second));
}
break;
}
case 3 : {
std::map<const std::string, std::set<MFace, Less_Face> >::iterator mapIt;
for(mapIt=boundaryFaces.begin(); mapIt!=boundaryFaces.end(); mapIt++) {
- bGroups.push_back(new dgGroupOfFaces(*eGroups[0],mapIt->first,order,mapIt->second));
+ bGroups.push_back(new dgGroupOfFaces(*eGroups[i],mapIt->first,order,mapIt->second));
}
break;
}
+ }
+ // create a group of faces for every face that is common to elements of the same group
+ fGroups.push_back(new dgGroupOfFaces(*eGroups[i],order));
+ // create a groupe of d
+ for (int j=i+1;j<eGroups.size();j++){
+ dgGroupOfFaces *gof = new dgGroupOfFaces(*eGroups[i],*eGroups[j],order);
+ if (gof->getNbElements())
+ fGroups.push_back(gof);
+ else
+ delete gof;
+ }
}
- fGroups.push_back(new dgGroupOfFaces(*eGroups[0],order));
}
-
// works for any number of groups
void dgAlgorithm::residual( const dgConservationLaw &claw,
std::vector<dgGroupOfElements*> &eGroups, //group of elements
@@ -616,3 +647,51 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
}
// residu[0]->print("Boundaries");
}
+
+void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF manager (maybe useless here)
+ const dgConservationLaw &claw, // the conservation law
+ const dgGroupOfElements & group,
+ fullMatrix<double> &solution,
+ std::vector<double> & DT
+ )
+{
+ dataCacheMap cacheMap(group.getNbNodes());
+ dataCacheDouble &sol = cacheMap.provideData("Solution");
+ dataCacheElement &cacheElement = cacheMap.getElement();
+ // provided dataCache
+ dataCacheDouble *maxConvectiveSpeed = claw.newMaxConvectiveSpeed(cacheMap);
+ dataCacheDouble *maximumDiffusivity = claw.newMaximumDiffusivity(cacheMap);
+
+ const int nbFields = claw.nbFields();
+ /* This is an estimate on how lengths changes with p
+ It is merely the smallest distance between gauss
+ points at order p + 1 */
+ const double p = group.getOrder();
+ const double Cip = 3 * (p + 1) * (p + group.getDimUVW()) ;
+ double l_red = 1./3.*p*p +7./6.*p +1.0;
+ double l_red_sq = l_red*l_red;
+ DT.resize(group.getNbElements());
+ for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
+ sol.setAsProxy(solution, iElement*nbFields, nbFields);
+ MElement *e = group.getElement(iElement);
+ cacheElement.set(e);
+ const double L = e->minEdge();
+ double spectralRadius = 0.0;
+ if (maximumDiffusivity){
+ double nu = (*maximumDiffusivity)(0,0);
+ for (int k=1;k<group.getNbNodes();k++) nu = std::max((*maximumDiffusivity)(k,0), nu);
+ spectralRadius += group.getDimUVW()*nu/(L*L)*std::max(l_red_sq,6.*l_red*Cip);
+ spectralRadius += 4.0*nu*l_red*Cip/(L*L);
+ }
+ if (maxConvectiveSpeed){
+ double c = (*maxConvectiveSpeed)(0,0);
+ for (int k=1;k<group.getNbNodes();k++) c = std::max((*maxConvectiveSpeed)(k,0), c);
+ spectralRadius += 4.0*c*l_red/L;
+ // printf("coucou %g %g %g %g\n",c,l_red,L,spectralRadius);
+ }
+ DT[iElement] = 1./spectralRadius;
+ }
+ delete maximumDiffusivity;
+ delete maxConvectiveSpeed;
+}
+
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index fe3bd7f36ad0ac60aa1eeeeb429aa192d24f5553..2402b2a6093ee6ea6889fb91a7b1531faae5ddc7 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -50,6 +50,12 @@ class dgAlgorithm {
dgDofContainer &sol,
dgLimiter *limiter=NULL,
int orderRK=4);
+ static void computeElementaryTimeSteps ( //dofManager &dof, // the DOF manager (maybe useless here)
+ const dgConservationLaw &claw, // the conservation law
+ const dgGroupOfElements & group, // the group
+ fullMatrix<double> &solution, // the solution
+ std::vector<double> & DT // elementary time steps
+ );
void multirateRungeKutta (
const dgConservationLaw &claw,
std::vector<dgGroupOfElements*> &eGroups, //group of elements
@@ -59,6 +65,7 @@ class dgAlgorithm {
dgDofContainer &residu,
dgDofContainer &sol,
int orderRK=3);
+
static void multAddInverseMassMatrix ( /*dofManager &dof,*/
const dgGroupOfElements & group,
fullMatrix<double> &residu,
diff --git a/Solver/dgConservationLawWaveEquation.cpp b/Solver/dgConservationLawWaveEquation.cpp
index 440c8993f7fb89d8770c40d45f914cd5e78db5a5..9d9e18ef7247aa99aa4fbb6dad67b8039760375c 100644
--- a/Solver/dgConservationLawWaveEquation.cpp
+++ b/Solver/dgConservationLawWaveEquation.cpp
@@ -28,6 +28,22 @@ class dgConservationLawWaveEquation::hyperbolicFlux : public dataCacheDouble {
}
}
};
+
+class dgConservationLawWaveEquation::maxConvectiveSpeed : public dataCacheDouble {
+ dataCacheDouble /
+ public:
+ maxConvectiveSpeed(dataCacheMap &cacheMap):
+ sol(cacheMap.get("Solution",this))
+ {
+ };
+ void _eval () {
+ int nQP = sol().size1();
+ if(_value.size1() != nQP)
+ _value=fullMatrix<double>(nQP,1);
+ _value.setAll(c);
+ }
+};
+
class dgConservationLawWaveEquation::riemann : public dataCacheDouble {
dataCacheDouble &normals, &solL, &solR;
const int _DIM,_nbf;
@@ -69,6 +85,9 @@ class dgConservationLawWaveEquation::riemann : public dataCacheDouble {
dataCacheDouble *dgConservationLawWaveEquation::newConvectiveFlux( dataCacheMap &cacheMap) const {
return new hyperbolicFlux(cacheMap,_DIM);
}
+dataCacheDouble *dgConservationLawWaveEquation::newMaxConvectiveSpeed( dataCacheMap &cacheMap) const {
+ return new maxConvectiveSpeed(cacheMap);
+}
dataCacheDouble *dgConservationLawWaveEquation::newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {
return new riemann(cacheMapLeft, cacheMapRight,_DIM);
}
diff --git a/Solver/dgConservationLawWaveEquation.h b/Solver/dgConservationLawWaveEquation.h
index adb3c2d43e94ce8de1d0162c25749bdabb71f880..530dff8b2d080b4f57c9e57f766fa4e06e16ded4 100644
--- a/Solver/dgConservationLawWaveEquation.h
+++ b/Solver/dgConservationLawWaveEquation.h
@@ -6,12 +6,14 @@ class constructorBinding;
class dgConservationLawWaveEquation : public dgConservationLaw {
int _DIM;
class hyperbolicFlux;
+ class maxConvectiveSpeed;
class riemann;
public:
dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const ;
dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const;
dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const ;
dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const ;
+ dataCacheDouble *newMaxConvectiveSpeed (dataCacheMap &cacheMap) const;
dgBoundaryCondition *newBoundaryWall () const;
dgConservationLawWaveEquation(int DIM);
static const char *className;
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 182e6652cb155947fe83f441f4529147a46a3323..c57c47a8ccb326325ef41b300b0aded3c9e8a02b 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -3,6 +3,7 @@
#include "polynomialBasis.h"
#include "Numeric.h"
#include "MTriangle.h"
+#include "MQuadrangle.h"
#include "MLine.h"
#include "MPoint.h"
#include "GModel.h"
@@ -153,11 +154,21 @@ void dgGroupOfFaces::addFace(const MFace &topoFace, int iElLeft, int iElRight){
// triangular face
if (topoFace.getNumVertices() == 3){
switch(vertices.size()){
- case 3 : _faces.push_back(new MTriangle (vertices) ); break;
- case 6 : _faces.push_back(new MTriangle6 (vertices) ); break;
- case 10 : _faces.push_back(new MTriangleN (vertices, 3) ); break;
- case 15 : _faces.push_back(new MTriangleN (vertices, 4) ); break;
- case 21 : _faces.push_back(new MTriangleN (vertices, 5) ); break;
+ case 3 : _faces.push_back(new MTriangle (vertices) ); break;
+ case 6 : _faces.push_back(new MTriangle6 (vertices) ); break;
+ case 10 : _faces.push_back(new MTriangleN (vertices,3) ); break;
+ case 15 : _faces.push_back(new MTriangleN (vertices,4) ); break;
+ case 21 : _faces.push_back(new MTriangleN (vertices,5) ); break;
+ default : throw;
+ }
+ }
+ else if (topoFace.getNumVertices() == 4){
+ switch(vertices.size()){
+ case 4 : _faces.push_back(new MQuadrangle (vertices) ); break;
+ case 8 : _faces.push_back(new MQuadrangle8 (vertices) ); break;
+ case 9 : _faces.push_back(new MQuadrangle9 (vertices) ); break;
+ case 16 : _faces.push_back(new MQuadrangleN (vertices, 4)); break;
+ case 25 : _faces.push_back(new MQuadrangleN (vertices, 5)); break;
default : throw;
}
}
@@ -182,11 +193,12 @@ void dgGroupOfFaces::addEdge(const MEdge &topoEdge, int iElLeft, int iElRight) {
const std::vector<int> &geomClosure = elLeft.getFunctionSpace()->getEdgeClosure(_closuresIdLeft.back());
for (int j=0; j<geomClosure.size() ; j++)
vertices.push_back( elLeft.getVertex(geomClosure[j]) );
- switch(vertices.size()){
+ switch(vertices.size())
+ {
case 2 : _faces.push_back(new MLine (vertices) ); break;
case 3 : _faces.push_back(new MLine3 (vertices) ); break;
default : _faces.push_back(new MLineN (vertices) ); break;
- }
+ }
}
void dgGroupOfFaces::addVertex(MVertex *topoVertex, int iElLeft, int iElRight){
@@ -205,6 +217,7 @@ void dgGroupOfFaces::addVertex(MVertex *topoVertex, int iElLeft, int iElRight){
}
void dgGroupOfFaces::init(int pOrder) {
+ if (!_faces.size())return;
_fsFace = _faces[0]->getFunctionSpace (pOrder);
_integration = dgGetIntegrationRule (_faces[0],pOrder);
_redistribution = new fullMatrix<double> (_fsFace->coefficients.size1(),_integration->size1());
@@ -299,6 +312,7 @@ void dgGroupOfFaces::init(int pOrder) {
dgGroupOfFaces::~dgGroupOfFaces()
{
+ if (!_faces.size())return;
delete _redistribution;
delete _collocation;
delete _detJac;
@@ -430,6 +444,107 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
init(pOrder);
}
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup1, const dgGroupOfElements &elGroup2, int pOrder):
+ _groupLeft(elGroup1),_groupRight(elGroup2)
+{
+ _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
+ _fsRight=_groupRight.getElement(0)->getFunctionSpace(pOrder);
+ switch (_groupLeft.getElement(0)->getDim()) {
+ case 1 : {
+ std::map<MVertex*,int> vertexMap1;
+ _closuresLeft = _fsLeft->vertexClosure;
+ _closuresRight = _fsRight->vertexClosure;
+
+ for(int i=0; i<elGroup1.getNbElements(); i++){
+ MElement &el = *elGroup1.getElement(i);
+ for (int j=0; j<el.getNumVertices(); j++){
+ MVertex* vertex = el.getVertex(j);
+ if(vertexMap1.find(vertex) == vertexMap1.end()){
+ vertexMap1[vertex] = i;
+ }else{
+ vertexMap1.erase(vertex);
+ }
+ }
+ }
+
+ for(int i=0; i<elGroup2.getNbElements(); i++){
+ MElement &el = *elGroup2.getElement(i);
+ for (int j=0; j<el.getNumVertices(); j++){
+ MVertex* vertex = el.getVertex(j);
+ std::map<MVertex*,int>::iterator it = vertexMap1.find(vertex);
+ if(it != vertexMap1.end()){
+ addVertex(vertex,it->second,i);
+ }
+ }
+ }
+ }
+ break;
+ case 2 : {
+ std::map<MEdge,int,Less_Edge> edgeMap;
+ _closuresLeft = _fsLeft->edgeClosure;
+ _closuresRight = _fsRight->edgeClosure;
+ for(int i=0; i<elGroup1.getNbElements(); i++){
+ MElement &el = *elGroup1.getElement(i);
+ for (int j=0; j<el.getNumEdges(); j++){
+ MEdge edge = el.getEdge(j);
+ if(edgeMap.find(edge) == edgeMap.end()){
+ edgeMap[edge] = i;
+ }else{
+ edgeMap.erase(edge);
+ }
+ }
+ }
+ for(int i=0; i<elGroup2.getNbElements(); i++){
+ MElement &el = *elGroup2.getElement(i);
+ for (int j=0; j<el.getNumEdges(); j++){
+ MEdge edge = el.getEdge(j);
+ std::map<MEdge,int,Less_Edge>::iterator it = edgeMap.find(edge);
+ if(it != edgeMap.end()){
+// MElement &el2 = *elGroup1.getElement(it->second);
+// printf("%d %d \n",edge.getVertex(0)->getNum(),edge.getVertex(1)->getNum());
+// for (int k=0;k<el.getNumVertices();k++)printf("%d ",el.getVertex(k)->getNum());
+// printf("\n");
+// for (int k=0;k<el2.getNumVertices();k++)printf("%d ",el2.getVertex(k)->getNum());
+// printf("\n");
+ addEdge(edge,it->second,i);
+ }
+ }
+ }
+
+ printf("groupe mixte %d %d elements %d faces\n",elGroup1.getNbElements(),elGroup2.getNbElements(),_faces.size());
+ break;
+ }
+ case 3 : {
+ std::map<MFace,int,Less_Face> faceMap;
+ _closuresLeft = _fsLeft->faceClosure;
+ _closuresRight = _fsRight->faceClosure;
+ for(int i=0; i<elGroup1.getNbElements(); i++){
+ MElement &el = *elGroup1.getElement(i);
+ for (int j=0; j<el.getNumFaces(); j++){
+ MFace face = el.getFace(j);
+ if(faceMap.find(face) == faceMap.end()){
+ faceMap[face] = i;
+ }else{
+ faceMap.erase(face);
+ }
+ }
+ }
+ for(int i=0; i<elGroup2.getNbElements(); i++){
+ MElement &el = *elGroup2.getElement(i);
+ for (int j=0; j<el.getNumFaces(); j++){
+ MFace face = el.getFace(j);
+ std::map<MFace,int,Less_Face>::iterator it = faceMap.find(face);
+ if(it != faceMap.end()){
+ addFace(face,it->second,i);
+ }
+ }
+ }
+ break;
+ }
+ default : throw;
+ }
+ init(pOrder);
+}
void dgGroupOfFaces::mapToInterface ( int nFields,
const fullMatrix<double> &vLeft,
const fullMatrix<double> &vRight,
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 224e71983fb9dae32c5d50e1b138648786b79cda..7d97db7081fbb150fe7103d645a6aee34a48382e 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -155,6 +155,7 @@ public:
inline fullMatrix<double> &getNormals () const {return *_normals;}
dgGroupOfFaces (const dgGroupOfElements &elements,int pOrder);
+ dgGroupOfFaces (const dgGroupOfElements &a, const dgGroupOfElements &b,int pOrder);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MVertex*> &boundaryVertices);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MEdge,Less_Edge> &boundaryEdges);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MFace,Less_Face> &boundaryFaces);
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index cb394b03402ff76cb3f877fa9f165afb41d31edd..b514d7a0e882e443be6a85b623e185d35b7ba662 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -50,9 +50,11 @@ methodBinding *dgSystemOfEquations::methods[]={
new methodBindingTemplate<dgSystemOfEquations,void>("limitSolution",&dgSystemOfEquations::limitSolution),
new methodBindingTemplate<dgSystemOfEquations,void,std::string>("L2Projection",&dgSystemOfEquations::L2Projection),
new methodBindingTemplate<dgSystemOfEquations,double,double>("RK44",&dgSystemOfEquations::RK44),
+ new methodBindingTemplate<dgSystemOfEquations,double>("computeInvSpectralRadius",&dgSystemOfEquations::computeInvSpectralRadius),
new methodBindingTemplate<dgSystemOfEquations,double,double>("RK44_limiter",&dgSystemOfEquations::RK44_limiter),
new methodBindingTemplate<dgSystemOfEquations,double,double>("multirateRK43",&dgSystemOfEquations::multirateRK43),
0};
+#endif // HAVE_LUA
// do a L2 projection
void dgSystemOfEquations::L2Projection (std::string functionName){
@@ -83,13 +85,23 @@ double dgSystemOfEquations::RK44(double dt){
return _solution->_data->norm();
}
+double dgSystemOfEquations::computeInvSpectralRadius(){
+
+ double sr = 1.e22;
+ for (int i=0;i<_elementGroups.size();i++){
+ std::vector<double> DTS;
+ _algo->computeElementaryTimeSteps(*_claw, *_elementGroups[i], *_solution->_dataProxys[i], DTS);
+ for (int k=0;k<DTS.size();k++) sr = std::min(sr,DTS[k]);
+ }
+ return sr;
+}
+
double dgSystemOfEquations::RK44_limiter(double dt){
dgLimiter *sl = new dgSlopeLimiter();
_algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide, sl);
return _solution->_data->norm();
}
-
double dgSystemOfEquations::multirateRK43(double dt){
_algo->multirateRungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide);
return _solution->_data->norm();
@@ -98,13 +110,13 @@ double dgSystemOfEquations::multirateRK43(double dt){
void dgSystemOfEquations::exportSolution(std::string outputFile){
export_solution_as_is(outputFile);
}
+
void dgSystemOfEquations::limitSolution(){
dgLimiter *sl = new dgSlopeLimiter();
sl->apply(*_solution,_elementGroups,_faceGroups);
delete sl;
}
-#endif // HAVE_LUA
dgSystemOfEquations::~dgSystemOfEquations(){
for (int i=0;i<_elementGroups.size();i++){
@@ -196,7 +208,7 @@ dgDofContainer::dgDofContainer (std::vector<dgGroupOfElements*> &elementGroups,
for (int i=0;i<elementGroups.size();i++){
int nbNodes = elementGroups[i]->getNbNodes();
int nbElements = elementGroups[i]->getNbElements();
- _dataProxys[i] = new fullMatrix<double> (&(*_data)(offset*sizeof(double)),nbNodes, nbFields*nbElements);
+ _dataProxys[i] = new fullMatrix<double> (&(*_data)(offset),nbNodes, nbFields*nbElements);
offset += nbNodes*nbFields*nbElements;
}
// printf("datasize = %d\n",_dataSize);
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index cdcb87a1eb13c96600d9d15cd741785ba275b162..32aa5d8c2e57f3ae292e8324f936d16894f55001 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -51,6 +51,7 @@ class dgSystemOfEquations {
// groups of faces (boundary conditions)
std::vector<dgGroupOfFaces*> _boundaryGroups;
dgSystemOfEquations(const dgSystemOfEquations &) {}
+ double computeTimeStepMethodOfLines () const;
public:
void setOrder (int order); // set the polynomial order
void setConservationLaw (dgConservationLaw *law); // set the conservationLaw
@@ -62,6 +63,7 @@ public:
double RK44_limiter (double dt);
double multirateRK43 (double dt);
void L2Projection (std::string functionName); // assign the solution to a given function
+ double computeInvSpectralRadius();
static const char className[];
static const char parentClassName[];