Skip to content
Snippets Groups Projects
Select Git revision
  • 04a6e4517080f516e787cb4b7fe2eda8aacb39db
  • 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

gmshProjection.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    FiniteElement.cpp 6.75 KiB
    // Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file for license information. Please report all
    // bugs and problems to <gmsh@geuz.org>.
    
    #include <stdlib.h>
    #include "GmshConfig.h"
    #include "GModel.h"
    #include "FiniteElement.h"
    #include "gmshAssembler.h"
    #include "gmshProjection.h"
    #include "gmshLaplace.h"
    #include "gmshHelmholtz.h"
    #include "gmshLinearSystemGmm.h"
    #include "mathEvaluator.h"
    
    StringXNumber FiniteElementOptions_Number[] = {
      {GMSH_FULLRC, "Omega", NULL, 1.},
      {GMSH_FULLRC, "Gamma1", NULL, 0.},
      {GMSH_FULLRC, "Gamma1Value", NULL, 0.},
      {GMSH_FULLRC, "Gamma2", NULL, 0.},
      {GMSH_FULLRC, "Gamma2Value", NULL, 0.}
    };
    
    StringXString FiniteElementOptions_String[] = {
      {GMSH_FULLRC, "Equation", NULL, "Projection"},
      {GMSH_FULLRC, "Parameter", NULL, "Sin(x*y)"},
      {GMSH_FULLRC, "Gamma1BC", NULL, ""},
      {GMSH_FULLRC, "Gamma2BC", NULL, ""},
    };
    
    extern "C"
    {
      GMSH_Plugin *GMSH_RegisterFiniteElementPlugin()
      {
        return new GMSH_FiniteElementPlugin();
      }
    }
    
    std::string GMSH_FiniteElementPlugin::getHelp() const
    {
      return "Plugin(FiniteElement) solves simple PDEs\n"
             "using the finite element method. This is only\n"
             "intended as a demonstration tool: it is NOT\n"
             "intended for general use."
             "\n"
             "Plugin(FiniteElement) creates a new view.\n";
    }
    
    int GMSH_FiniteElementPlugin::getNbOptions() const
    {
      return sizeof(FiniteElementOptions_Number) / sizeof(StringXNumber);
    }
    
    StringXNumber *GMSH_FiniteElementPlugin::getOption(int iopt)
    {
      return &FiniteElementOptions_Number[iopt];
    }
    
    int GMSH_FiniteElementPlugin::getNbOptionsStr() const
    {
      return sizeof(FiniteElementOptions_String) / sizeof(StringXString);
    }
    
    StringXString *GMSH_FiniteElementPlugin::getOptionStr(int iopt)
    {
      return &FiniteElementOptions_String[iopt];
    }
    
    template<class scalar>
    class gmshMathEvalFunction : public simpleFunction<scalar> {
     private:
      std::string _str;
      mathEvaluator *_f;
     public:
      gmshMathEvalFunction(std::string str) : _str(str)
      {
        std::vector<std::string> expressions(1), variables(3);
        expressions[0] = str;
        variables[0] = "x";
        variables[1] = "y";
        variables[2] = "z";
        _f = new mathEvaluator(expressions, variables);
        if(expressions.empty()){
          delete _f;
          _f = 0;
        }
      }
      virtual ~gmshMathEvalFunction()
      {
        if(_f) delete _f;
      }
      virtual scalar operator () (double x, double y, double z) const 
      { 
        if(_f){
          std::vector<double> values(3), res(1);
          values[0] = x;
          values[1] = y;
          values[2] = z;
          if(_f->eval(values, res)) return res[0];
        }
        return atof(_str.c_str());
      }
    };
    
    #if defined(HAVE_GMM)
    template<class scalar>
    class solver{
     public:
      gmshLinearSystemGmm<scalar> *lsys;
      gmshAssembler<scalar> *myAssembler;
      solver()
      {
        int volume = (int)FiniteElementOptions_Number[0].def;
        int surface1 = (int)FiniteElementOptions_Number[1].def;
        double value1 = FiniteElementOptions_Number[2].def;
        int surface2 = (int)FiniteElementOptions_Number[3].def;
        double value2 = FiniteElementOptions_Number[4].def;
        std::string bc1 = FiniteElementOptions_String[2].def;
        std::string bc2 = FiniteElementOptions_String[3].def;
    
        lsys = new gmshLinearSystemGmm<scalar>;
        lsys->setPrec(1.e-10);
        lsys->setNoisy(2);
        myAssembler = new gmshAssembler<scalar>(lsys);
    
        GModel *m = GModel::current();
        std::vector<MVertex*> vertices;
    
        int dim = m->getNumRegions() ? 3 : 2;
        
        if(bc1 == "Dirichlet"){
          m->getMeshVerticesForPhysicalGroup(dim - 1, surface1, vertices);
          for(unsigned int i = 0; i < vertices.size(); i++)
            myAssembler->fixVertex(vertices[i], 0, 1, value1);
        }
        if(bc2 == "Dirichlet"){
          m->getMeshVerticesForPhysicalGroup(dim - 1, surface2, vertices);
          for(unsigned int i = 0; i < vertices.size(); i++)
            myAssembler->fixVertex(vertices[i], 0, 1, value2);
        }
        m->getMeshVerticesForPhysicalGroup(dim, volume, vertices);
        for(unsigned int i = 0; i < vertices.size(); i++)
          myAssembler->numberVertex(vertices[i], 0, 1);
    
        Msg::Info("Assembler: %d dofs, %d fixed",
                  myAssembler->sizeOfR(), myAssembler->sizeOfF());
      }
    };
    #endif
    
    PView *GMSH_FiniteElementPlugin::execute(PView *v)
    {
      std::string equation = FiniteElementOptions_String[0].def;
      std::string parameter = FiniteElementOptions_String[1].def;
      int volume = (int)FiniteElementOptions_Number[0].def;
    
    #if defined(HAVE_GMM)
      GModel *m = GModel::current();
      std::map<int, std::vector<GEntity*> > groups[4];
      m->getPhysicalGroups(groups);
      int dim = m->getNumRegions() ? 3 : 2;
    
      std::vector<MVertex*> vertices;
      m->getMeshVerticesForPhysicalGroup(dim, volume, vertices);
      std::map<int, std::vector<double> > data;
    
      if(equation == "Projection"){
        solver<double> s;
        if(!s.myAssembler->sizeOfR()) return 0;
        gmshMathEvalFunction<double> par(parameter);
        gmshProjectionTerm projection(m, 1);
        for(unsigned int i = 0; i < groups[dim][volume].size(); i++)
          projection.addToMatrix(*s.myAssembler, groups[dim][volume][i]);
        projection.addNeumann(volume, dim, 0, 1, par, *s.myAssembler);
        s.lsys->systemSolve();
        for (unsigned int i = 0; i < vertices.size(); i++)
          data[vertices[i]->getNum()].push_back
            (s.myAssembler->getDofValue(vertices[i], 0, 1));
        new PView("projection", "NodeData", m, data, 0.);
      }
      else if(equation == "Laplace"){
        solver<double> s;
        if(!s.myAssembler->sizeOfR()) return 0;
        gmshMathEvalFunction<double> par(parameter);
        gmshLaplaceTerm laplace(m, &par, 1);
        for(unsigned int i = 0; i < groups[dim][volume].size(); i++)
          laplace.addToMatrix(*s.myAssembler, groups[dim][volume][i]);
        s.lsys->systemSolve();
        for (unsigned int i = 0; i < vertices.size(); i++)
          data[vertices[i]->getNum()].push_back
            (s.myAssembler->getDofValue(vertices[i], 0, 1));
        new PView("laplace", "NodeData", m, data, 0.);
      }
      else if(equation == "Helmholtz"){
        solver<std::complex<double> > s;
        if(!s.myAssembler->sizeOfR()) return 0;
        gmshMathEvalFunction<std::complex<double> > par(parameter);
        gmshHelmholtzTerm helmholtz(m, &par, 1);
        for(unsigned int i = 0; i < groups[dim][volume].size(); i++)
          helmholtz.addToMatrix(*s.myAssembler, groups[dim][volume][i]);
        s.lsys->setGmres(1);
        s.lsys->systemSolve();
        for (unsigned int i = 0; i < vertices.size(); i++)
          data[vertices[i]->getNum()].push_back
            (s.myAssembler->getDofValue(vertices[i], 0, 1).real());
        PView *pv = new PView("helmholtz", "NodeData", m, data, 0.);
        data.clear();
        for (unsigned int i = 0; i < vertices.size(); i++)
          data[vertices[i]->getNum()].push_back
            (s.myAssembler->getDofValue(vertices[i], 0, 1).imag());
        pv->addStep(m, data, 1.);
      }
      else{
        Msg::Error("Unknown equation to solve");
      }
    #endif
      return 0;
    }