From 87bdce186428110875bcfb8d3ea6320a174df5f0 Mon Sep 17 00:00:00 2001
From: Samuel Melchior <samuel.melchior@uclouvain.be>
Date: Tue, 8 Dec 2009 08:00:39 +0000
Subject: [PATCH] dgMainSam
---
Solver/dgMainSam.cpp | 110 +++++++++++++++++++++++++++++++++++++++++++
1 file changed, 110 insertions(+)
create mode 100644 Solver/dgMainSam.cpp
diff --git a/Solver/dgMainSam.cpp b/Solver/dgMainSam.cpp
new file mode 100644
index 0000000000..a33a33331a
--- /dev/null
+++ b/Solver/dgMainSam.cpp
@@ -0,0 +1,110 @@
+#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);
+}
--
GitLab