Skip to content
Snippets Groups Projects
Select Git revision
  • 07fe1c9c3ec4bae605653baf3890d51c5b047115
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
  • gmsh_2_8_6
26 results

PViewDataGModel.h

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    dgMainTest.cc 3.57 KiB
    #include <stdio.h>
    #include <vector>
    #include "GModel.h"
    #include "dgGroupOfElements.h"
    #include "dgAlgorithm.h"
    #include "dgConservationLaw.h"
    #include "Gmsh.h"
    
    
    #include "MElement.h"
    void print (const char *filename,const dgGroupOfElements &els, double *v);
    std::vector<MElement *> getAllTri(GModel *model);
    
    class testSourceTerm : public dgTerm {
      void operator () (const dgElement &el, fullMatrix<double> fcx[]) const{
        const fullMatrix<double> &sol = el.solution();
        const fullMatrix<double> &qp = el.integration();
        SPoint3 p;
        for(int i=0; i< sol.size1(); i++) {
          el.element()->pnt(qp(i,0),qp(i,1),qp(i,2),p);
          fcx[0](i,0)=exp(-(pow(p.x()-0.2,2)+pow(p.y(),2))*100);
        }
      }
    };
    
    class dgConservationLawInitialCondition : public dgConservationLaw {
      public:
      dgConservationLawInitialCondition() {
        _nbf = 1;
        _source = new testSourceTerm;
      }
      ~dgConservationLawInitialCondition() {
        delete _source;
      }
    };
    
    int main(int argc, char **argv){
      GmshMergeFile("input/mesh1.msh");
      //GmshMergeFile("square2.msh");
      std::vector<MElement *> allTri=getAllTri(GModel::current());
      int order=1;
      dgGroupOfElements elements(allTri,order);
      dgGroupOfFaces faces(elements,order);
      //dgGroupOfFaces faces(elements0,elements1,order);
      int nbNodes = elements.getNbNodes();
      fullMatrix<double> sol(nbNodes,elements.getNbElements());
      fullMatrix<double> residu(nbNodes,elements.getNbElements());
      dgAlgorithm algo;
      // initial condition
      {
        dgConservationLawInitialCondition initLaw;
        algo.residualVolume(initLaw,elements,sol,residu);
        algo.multAddInverseMassMatrix(elements,residu,sol);
      }
      print("init.pos",elements,&sol(0,0));
      //advection
      dgConservationLaw *law = dgNewConservationLawAdvection();
      for(int iT=0; iT<1000; iT++) {
        residu.scale(0);
        algo.residualVolume(*law,elements,sol,residu);
        { //interface term
          fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2);
          fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2);
          faces.mapToInterface(1, sol, sol, solInterface);
          algo.residualInterface(*law,faces,solInterface,residuInterface);
          faces.mapFromInterface(1, residuInterface, residu, residu);
        }
        residu.scale(0.01);
        algo.multAddInverseMassMatrix(elements,residu,sol);
        if(iT%10==0){
          char name[100];
          sprintf(name,"test_%05i.pos",iT/10);
          printf("%i\n",iT);
          print(name,elements,&sol(0,0));
        }
      }
      delete law;
    }
    
    std::vector<MElement *> getAllTri(GModel *model){
      std::vector<GEntity*> entities;
      model->getEntities(entities);
      std::vector<MElement *> all_tri;
      for(std::vector<GEntity *>::iterator itent=entities.begin(); itent!=entities.end(); itent++){
        if ((*itent)->dim()!=2) continue;
        for (int iel=0; iel<(*itent)->getNumMeshElements(); iel++)
          all_tri.push_back((*itent)->getMeshElement(iel));
      }
      return all_tri;
    }
    
    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,"ST (");
        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);
    }