diff --git a/CMakeLists.txt b/CMakeLists.txt index dba191531708f11a42bc2b2dc2d80a862cef508e..0e30ba2a9a0365481a9b81ee578f86d4ea9ccd0d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1007,6 +1007,3 @@ mark_as_advanced(BISON FLEX GMP_LIB GMSH_EXTRA_VERSION HDF5_LIB MAKEINFO add_executable(gmshlua EXCLUDE_FROM_ALL Solver/dgMainLua.cpp ${GMSH_SRC}) target_link_libraries(gmshlua ${LINK_LIBRARIES}) - -add_executable(gmshdgsam EXCLUDE_FROM_ALL Solver/dgMainSam.cpp ${GMSH_SRC}) -target_link_libraries(gmshdgsam ${LINK_LIBRARIES}) diff --git a/Solver/dgMainSam.cpp b/Solver/dgMainSam.cpp deleted file mode 100644 index a33a33331a66ff6e6b088f3180f97757001c9060..0000000000000000000000000000000000000000 --- a/Solver/dgMainSam.cpp +++ /dev/null @@ -1,110 +0,0 @@ -#include <stdio.h> -#include <vector> -#include "GModel.h" -#include "dgGroupOfElements.h" -#include "dgAlgorithm.h" -#include "dgConservationLaw.h" -#include "Gmsh.h" -#include "function.h" - - -#include "MElement.h" -void print (const char *filename,const dgGroupOfElements &els, double *v); - - -class dgConservationLawInitialCondition : public dgConservationLaw { - class gaussian : public dataCacheDouble { - dataCacheDouble &xyz; - double _xc; - public: - gaussian(dataCacheMap &cacheMap,double xc):xyz(cacheMap.get("XYZ",this)),_xc(xc){}; - void _eval () { - if(_value.size1() != xyz().size1()) - _value=fullMatrix<double>(xyz().size1(),1); - for(int i=0; i< _value.size1(); i++) { - _value(i,0)= exp(-pow(xyz(i,0)-_xc,2)*100); //1.0-xyz(i,0);// - } - } - }; - public: - dgConservationLawInitialCondition() { - _nbf = 1; - } - dataCacheDouble *newSourceTerm(dataCacheMap &cacheMap)const { - return new gaussian(cacheMap,-0.8); - } -}; - -int main(int argc, char **argv){ - GmshMergeFile("input/edge.msh"); - //we probably need a class to group those three ones - std::vector<dgGroupOfElements*> elementGroups; - std::vector<dgGroupOfFaces*> faceGroups; - std::vector<dgGroupOfFaces*> boundaryGroups; - int order=1; - int dimension=1; - dgAlgorithm algo; - function::registerDefaultFunctions(); - algo.buildGroups(GModel::current(),dimension,order,elementGroups,faceGroups,boundaryGroups); - - //for now, we suppose there is only one group of elements - int nbNodes = elementGroups[0]->getNbNodes(); - fullMatrix<double> sol(nbNodes,elementGroups[0]->getNbElements()); - fullMatrix<double> residu(nbNodes,elementGroups[0]->getNbElements()); - // initial condition - { - dgConservationLawInitialCondition initLaw; - algo.residualVolume(initLaw,*elementGroups[0],sol,residu); - algo.multAddInverseMassMatrix(*elementGroups[0],residu,sol); - } - - - - - fullMatrix<double> advectionSpeed(3,1); - advectionSpeed(0,0)=1.0; - advectionSpeed(1,0)=0; - advectionSpeed(2,0)=0.; - function::add("advectionSpeed",function::newFunctionConstant(advectionSpeed)); - - dgConservationLaw *law = dgNewConservationLawAdvection("advectionSpeed"); - - fullMatrix<double> outsideValue(1,1); - outsideValue(0,0)=0; - function::add("outsideValue",function::newFunctionConstant(outsideValue)); - law->addBoundaryCondition("Left",dgBoundaryCondition::newOutsideValueCondition(*law,"outsideValue")); - law->addBoundaryCondition("Right",dgBoundaryCondition::newOutsideValueCondition(*law,"outsideValue")); - - print("output/init.pos",*elementGroups[0],&sol(0,0)); - int n_visu = 10; - for(int iT=0; iT<100*n_visu; iT++) { - algo.rungeKutta(*law,elementGroups,faceGroups,boundaryGroups,0.002,residu,sol); - if(iT%n_visu==0){ - char name[100]; - sprintf(name,"output/test_%05i.pos",iT/n_visu); - printf("%i\n",iT); - print(name,*elementGroups[0],&sol(0,0)); - } - } - delete law; -} - -void print (const char *filename,const dgGroupOfElements &els, double *v) { - FILE *file = fopen(filename,"w"); - fprintf(file,"View \"%s\" {\n", filename); - for(int iel=0;iel<els.getNbElements();iel++){ - MElement *el = els.getElement(iel); - fprintf(file,"SL ("); - for (int iv=0; iv<el->getNumVertices(); iv++) { - MVertex *vertex = el->getVertex(iv); - SPoint3 coord = vertex->point(); - fprintf(file,"%e, %e, %e%c ",coord.x(),coord.y(),coord.z(),iv==el->getNumVertices()-1?')':','); - } - fprintf(file,"{"); - for (int iv=0; iv<el->getNumVertices(); iv++) - fprintf(file,"%e%c ",*(v++),iv==el->getNumVertices()-1?'}':','); - fprintf(file,";\n"); - } - fprintf(file,"};"); - fclose(file); -}