Skip to content
Snippets Groups Projects
Select Git revision
  • 3f147e49852e19b79b873f3f88edaff2314ae196
  • master default protected
  • overlaps_tags_and_distributed_export
  • overlaps_tags_and_distributed_export_rebased
  • relaying
  • alphashapes
  • patches-4.14
  • steplayer
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
41 results

solverAlgorithms.h

Blame
  • solverAlgorithms.h 10.53 KiB
    // Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file for license information. Please report all
    // issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
    //
    // Contributor(s):
    //   Eric Bechet
    //
    
    #ifndef SOLVERALGORITHMS_H
    #define SOLVERALGORITHMS_H
    
    #include "dofManager.h"
    #include "terms.h"
    #include "quadratureRules.h"
    #include "MVertex.h"
    
    template <class Iterator, class Assembler>
    void Assemble(BilinearTermBase &term, FunctionSpaceBase &space,
                  Iterator itbegin, Iterator itend, QuadratureBase &integrator,
                  Assembler &assembler)
    // symmetric
    {
      fullMatrix<typename Assembler::dataMat> localMatrix;
      std::vector<Dof> R;
      for(Iterator it = itbegin; it != itend; ++it) {
        MElement *e = *it;
        R.clear();
        IntPt *GP;
        int npts = integrator.getIntPoints(e, &GP);
        term.get(e, npts, GP, localMatrix); // localMatrix.print();
        space.getKeys(e, R);
        assembler.assemble(R, localMatrix);
      }
    }
    
    template <class Iterator, class Assembler>
    void Assemble(BilinearTermBase &term, FunctionSpaceBase &space,
                  Iterator itbegin, Iterator itend, QuadratureBase &integrator,
                  Assembler &assembler, elementFilter &efilter)
    // symmetric
    {
      fullMatrix<typename Assembler::dataMat> localMatrix;
      std::vector<Dof> R;
      for(Iterator it = itbegin; it != itend; ++it) {
        MElement *e = *it;
        if(efilter(e)) {
          R.clear();
          IntPt *GP;
          int npts = integrator.getIntPoints(e, &GP);
          term.get(e, npts, GP, localMatrix); // localMatrix.print();
          space.getKeys(e, R);
          assembler.assemble(R, localMatrix);
        }
      }
    }
    
    template <class Assembler>
    void Assemble(BilinearTermBase &term, FunctionSpaceBase &space, MElement *e,
                  QuadratureBase &integrator, Assembler &assembler) // symmetric
    {
      fullMatrix<typename Assembler::dataMat> localMatrix;
      std::vector<Dof> R;
      IntPt *GP;
      int npts = integrator.getIntPoints(e, &GP);
      term.get(e, npts, GP, localMatrix);
      space.getKeys(e, R);
      assembler.assemble(R, localMatrix);
    }
    
    template <class Iterator, class Assembler>
    void Assemble(BilinearTermBase &term, FunctionSpaceBase &shapeFcts,
                  FunctionSpaceBase &testFcts, Iterator itbegin, Iterator itend,
                  QuadratureBase &integrator,
                  Assembler &assembler) // non symmetric
    {
      fullMatrix<typename Assembler::dataMat> localMatrix;
      std::vector<Dof> R;
      std::vector<Dof> C;
      for(Iterator it = itbegin; it != itend; ++it) {
        MElement *e = *it;
        R.clear();
        C.clear();
        IntPt *GP;
        int npts = integrator.getIntPoints(e, &GP);
        term.get(e, npts, GP, localMatrix); // localMatrix.print();
        // printf("local matrix size = %d %d\n", localMatrix.size1(),
        // localMatrix.size2());
        shapeFcts.getKeys(e, R);
        testFcts.getKeys(e, C);
        //    std::cout << "assembling normal test function ; lagrange trial
        //    function : " << std::endl; for (int i = 0 ; i < R.size() ; i++)
        //    {
        //      std::cout << "tests : " << R[i].getEntity() << ":" << R[i].getType()
        //      << std::endl ;
        //    }
        //    for (int i = 0 ; i < R.size() ; i++)
        //    {
        //      std::cout << "trial : " << C[i].getEntity() << ":" << C[i].getType()
        //      << std::endl ;
        //    }
        assembler.assemble(R, C, localMatrix);
        //    std::cout << "assembling lagrange test function ; normal trial
        //    function : " << std::endl; for (int i = 0 ; i < R.size() ; i++)
        //    {
        //      std::cout << "tests : " << C[i].getEntity() << ":" << C[i].getType()
        //      << std::endl ;
        //    }
        //    for (int i = 0 ; i < R.size() ; i++)
        //    {
        //      std::cout << "trial : " << R[i].getEntity() << ":" << R[i].getType()
        //      << std::endl ;
        //    }
        assembler.assemble(C, R, localMatrix.transpose());
      }
    }
    
    template <class Iterator, class Assembler>
    void Assemble(LinearTermBase<double> &term, FunctionSpaceBase &space,
                  Iterator itbegin, Iterator itend, QuadratureBase &integrator,
                  Assembler &assembler)
    {
      fullVector<typename Assembler::dataMat> localVector;
      std::vector<Dof> R;
      for(Iterator it = itbegin; it != itend; ++it) {
        MElement *e = *it;
        R.clear();
        IntPt *GP;
        int npts = integrator.getIntPoints(e, &GP);
        term.get(e, npts, GP, localVector); // localVector.print();
        space.getKeys(e, R);
        assembler.assemble(R, localVector);
      }
    }
    
    template <class Iterator, class Assembler>
    void Assemble(LinearTermBase<double> &term, FunctionSpaceBase &space,
                  Iterator itbegin, Iterator itend, QuadratureBase &integrator,
                  Assembler &assembler, elementFilter &efilter)
    {
      fullVector<typename Assembler::dataMat> localVector;
      std::vector<Dof> R;
      for(Iterator it = itbegin; it != itend; ++it) {
        MElement *e = *it;
        if(efilter(e)) {
          R.clear();
          IntPt *GP;
          int npts = integrator.getIntPoints(e, &GP);
          term.get(e, npts, GP, localVector); // localVector.print();
          space.getKeys(e, R);
          assembler.assemble(R, localVector);
        }
      }
    }
    
    template <class Assembler>
    void Assemble(LinearTermBase<double> &term, FunctionSpaceBase &space,
                  MElement *e, QuadratureBase &integrator, Assembler &assembler)
    {
      fullVector<typename Assembler::dataMat> localVector;
      std::vector<Dof> R;
      IntPt *GP;
      int npts = integrator.getIntPoints(e, &GP);
      term.get(e, npts, GP, localVector);
      space.getKeys(e, R);
      assembler.assemble(R, localVector);
    }
    
    template <class Iterator, class dataMat>
    void Assemble(ScalarTermBase<double> &term, Iterator itbegin, Iterator itend,
                  QuadratureBase &integrator, dataMat &val)
    {
      dataMat localval;
      for(Iterator it = itbegin; it != itend; ++it) {
        MElement *e = *it;
        IntPt *GP;
        int npts = integrator.getIntPoints(e, &GP);
        term.get(e, npts, GP, localval);
        val += localval;
      }
    }
    
    template <class Iterator, class dataMat>
    void Assemble(ScalarTermBase<double> &term, Iterator itbegin, Iterator itend,
                  QuadratureBase &integrator, dataMat &val, elementFilter &efilter)
    {
      dataMat localval;
      for(Iterator it = itbegin; it != itend; ++it) {
        MElement *e = *it;
        if(efilter(e)) {
          IntPt *GP;
          int npts = integrator.getIntPoints(e, &GP);
          term.get(e, npts, GP, localval);
          val += localval;
        }
      }
    }
    
    template <class Iterator, class dataMat>
    void Assemble(ScalarTermBase<double> &term, MElement *e,
                  QuadratureBase &integrator, dataMat &val)
    {
      dataMat localval;
      IntPt *GP;
      int npts = integrator.getIntPoints(e, &GP);
      term.get(e, npts, GP, localval);
      val += localval;
    }
    
    template <class Assembler>
    void FixDofs(Assembler &assembler, std::vector<Dof> &dofs,
                 std::vector<typename Assembler::dataVec> &vals)
    {
      int nbff = dofs.size();
      for(int i = 0; i < nbff; ++i) {
        assembler.fixDof(dofs[i], vals[i]);
      }
    }
    
    class FilterDof {
    public:
      virtual ~FilterDof() {}
      virtual bool operator()(Dof key) = 0;
    };
    
    class FilterDofTrivial : public FilterDof {
    public:
      virtual bool operator()(Dof key) { return true; }
    };
    
    class FilterDofComponent : public FilterDof {
      int comp;
    
    public:
      FilterDofComponent(int comp_) : comp(comp_) {}
      virtual bool operator()(Dof key)
      {
        int type = key.getType();
        int icomp, iphys;
        Dof::getTwoIntsFromType(type, icomp, iphys);
        if(icomp == comp) return true;
        return false;
      }
    };
    
    // return true if the Dof is in a filled set
    class FilterDofSet : public FilterDof {
    protected:
      std::set<Dof> _dofset;
    
    public:
      FilterDofSet() : FilterDof() {}
      virtual ~FilterDofSet() {}
      virtual bool operator()(Dof key)
      {
        std::set<Dof>::iterator itR = _dofset.find(key);
        if(itR == _dofset.end()) {
          return false;
        }
        return true;
      }
      virtual void addDof(Dof key) { _dofset.insert(key); }
      virtual void addDof(std::vector<Dof> &R)
      {
        for(std::size_t i = 0; i < R.size(); i++) this->addDof(R[i]);
      }
    };
    
    template <class Assembler>
    void FixNodalDofs(FunctionSpaceBase &space, MElement *e, Assembler &assembler,
                      simpleFunction<typename Assembler::dataVec> &fct,
                      FilterDof &filter)
    {
      std::vector<MVertex *> tabV;
      int nv = e->getNumVertices();
      std::vector<Dof> R;
      space.getKeys(e, R);
      tabV.reserve(nv);
      for(int i = 0; i < nv; ++i) tabV.push_back(e->getVertex(i));
    
      for(std::vector<Dof>::iterator itd = R.begin(); itd != R.end(); ++itd) {
        Dof key = *itd;
        if(filter(key)) {
          for(int i = 0; i < nv; ++i) {
            if((long int)tabV[i]->getNum() == key.getEntity()) {
              assembler.fixDof(key, fct(tabV[i]->x(), tabV[i]->y(), tabV[i]->z()));
              break;
            }
          }
        }
      }
    }
    
    template <class Iterator, class Assembler>
    void FixNodalDofs(FunctionSpaceBase &space, Iterator itbegin, Iterator itend,
                      Assembler &assembler,
                      simpleFunction<typename Assembler::dataVec> &fct,
                      FilterDof &filter)
    {
      for(Iterator it = itbegin; it != itend; ++it)
        FixNodalDofs(space, *it, assembler, fct, filter);
    }
    
    template <class Iterator, class Assembler>
    void FixVoidNodalDofs(FunctionSpaceBase &space, Iterator itbegin,
                          Iterator itend, Assembler &assembler)
    {
      FilterDofTrivial filter;
      simpleFunction<typename Assembler::dataVec> fct(0.);
      FixNodalDofs(space, itbegin, itend, assembler, fct, filter);
    }
    
    template <class Iterator, class Assembler>
    void NumberDofs(FunctionSpaceBase &space, Iterator itbegin, Iterator itend,
                    Assembler &assembler)
    {
      for(Iterator it = itbegin; it != itend; ++it) {
        MElement *e = *it;
        std::vector<Dof> R;
        space.getKeys(e, R);
        int nbdofs = R.size();
        for(int i = 0; i < nbdofs; ++i) assembler.numberDof(R[i]);
      }
    }
    
      //// Mean HangingNodes
      // template <class Assembler> void FillHangingNodes(FunctionSpaceBase &space,
      // std::map<int,std::vector <int> > &HangingNodes, Assembler &assembler, int
      // &field, int &dim)
      //{
      //  std::map<int, std::vector <int> >::iterator ith;
      //  ith = HangingNodes.begin();
      //  int compt = 1;
      //  while (ith != HangingNodes.end()){
      //    float fac;
      //    fac = 1.0 / (ith->second).size();
      //    for (int j = 0; j < dim; j++){
      //      DofAffineConstraint<double> constraint;
      //      int type = Dof::createTypeWithTwoInts(j, field);
      //      Dof hgnd(ith->first, type);
      //      for (int i = 0; i < (ith->second).size(); i++){
      //          Dof key((ith->second)[i], type);
      //          std::pair<Dof, double > linDof(key, fac);
      //          constraint.linear.push_back(linDof);
      //      }
      //      constraint.shift = 0;
      //      assembler.setLinearConstraint (hgnd, constraint);
      //    }
      //    ith++;
      //    compt++;
      //  }
      //}
    
    #endif