Skip to content
Snippets Groups Projects
Select Git revision
  • 6d2114aac14613bceb69d488811c078870daf192
  • master default protected
  • hierarchical-basis
  • alphashapes
  • bl
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • convert_fdivs
  • tmp_jcjc24
  • fixedMeshIF
  • save_edges
  • 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

Numeric.cpp

Blame
  • frameSolver.cpp 9.24 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.
    
    #include "GmshConfig.h"
    #include "GModel.h"
    #include "GVertex.h"
    #include "GEdge.h"
    #include "frameSolver.h"
    #include "linearSystemCSR.h"
    #include "linearSystemPETSc.h"
    #include "linearSystemFull.h"
    
    #if defined(HAVE_POST)
    #include "PView.h"
    #include "PViewData.h"
    #endif
    
    frameSolver2d::frameSolver2d(GModel *gm) : pAssembler(0), _myModel(gm) {}
    
    void frameSolver2d::addFixations(const std::vector<int> &dirs,
                                     const std::vector<int> &modelVertices,
                                     double value)
    {
      for(std::size_t j = 0; j < modelVertices.size(); j++) {
        GVertex *gv = _myModel->getVertexByTag(modelVertices[j]);
        if(gv) {
          for(std::size_t i = 0; i < dirs.size(); i++) {
            _fixations.push_back(gmshFixation(gv, dirs[i], value));
          }
        }
      }
    }
    
    void frameSolver2d::addNodalForces(const std::vector<int> &modelVertices,
                                       const std::vector<double> &force)
    {
      for(std::size_t j = 0; j < modelVertices.size(); j++) {
        GVertex *gv = _myModel->getVertexByTag(modelVertices[j]);
        if(gv) {
          _nodalForces.push_back(std::make_pair(gv, force));
        }
      }
    }
    
    void frameSolver2d::addBeamsOrBars(const std::vector<int> &modelEdges, double E,
                                       double I, double A, int r[2])
    {
      int r_middle[2] = {1, 1}, r_left[2] = {r[0], 1}, r_right[2] = {0, r[1]};
      //  printf("adding %d beams\n",modelEdges.size());
      for(std::size_t i = 0; i < modelEdges.size(); i++) {
        GEdge *ge = _myModel->getEdgeByTag(modelEdges[i]);
        if(ge) {
          //      printf("model edge %d found\n",ge->tag());
          for(std::size_t j = 0; j < ge->lines.size(); ++j) {
            MLine *l = ge->lines[j];
            if(j == 0 && j == ge->lines.size() - 1)
              _beams.push_back(gmshBeam2d(l, E, I, A, r));
            else if(j == 0)
              _beams.push_back(gmshBeam2d(l, E, I, A, r_left));
            else if(j == ge->lines.size() - 1)
              _beams.push_back(gmshBeam2d(l, E, I, A, r_right));
            else
              _beams.push_back(gmshBeam2d(l, E, I, A, r_middle));
          }
        }
      }
    }
    
    void frameSolver2d::addBeams(const std::vector<int> &modelEdges, double E,
                                 double I, double A)
    {
      int r[2] = {1, 1};
      addBeamsOrBars(modelEdges, E, I, A, r);
    }
    
    void frameSolver2d::addBars(const std::vector<int> &modelEdges, double E,
                                double I, double A)
    {
      int r[2] = {0, 0};
      addBeamsOrBars(modelEdges, E, I, A, r);
    }
    
    // solve any time a parameter is modified
    /*
    
      A vertex is connected to beams. The question
      is how many bars are rotuled
    
      We define 2 dofs per node
    
     */
    void frameSolver2d::createDofs()
    {
      //  printf("coucou %d fixations\n",_fixations.size());
      for(std::size_t i = 0; i < _fixations.size(); ++i) {
        const gmshFixation &f = _fixations[i];
        //    printf("f._vertex(%d) = %p %d
        //    %g\n",i,f._vertex,f._direction,f._value);
        MVertex *v = f._vertex->mesh_vertices[0];
        Dof DOF(v->getNum(), f._direction);
        pAssembler->fixDof(DOF, f._value);
      }
    
      //  printf("coucou2\n");
      computeRotationTags();
      //  printf("coucou3\n");
      for(std::size_t i = 0; i < _beams.size(); i++) {
        //    printf("beam[%d] rot %d
        //    %d\n",i,_beams[i]._rotationTags[0],_beams[i]._rotationTags[1]);
        for(std::size_t j = 0; j < 2; j++) {
          MVertex *v = _beams[i]._element->getVertex(j);
          Dof theta(v->getNum(),
                    Dof::createTypeWithTwoInts(2, _beams[i]._rotationTags[j]));
          pAssembler->numberDof(theta);
          Dof U(v->getNum(), 0);
          pAssembler->numberDof(U);
          Dof V(v->getNum(), 1);
          pAssembler->numberDof(V);
        }
      }
      //  printf("%d dofs\n",pAssembler->sizeOfR());
    }
    
    void frameSolver2d::computeStiffnessMatrix(int iBeam, fullMatrix<double> &K)
    {
      const gmshBeam2d &b = _beams[iBeam];
      const double BS = b._e * b._i / (b._l * b._l * b._l);
      const double TS = b._e * b._a / b._l;
      const MVertex *v1 = b._element->getVertex(0);
      const MVertex *v2 = b._element->getVertex(1);
      const double alpha = atan2(v2->y() - v1->y(), v2->x() - v1->x());
      const double C = cos(alpha);
      const double S = sin(alpha);
    
      printf("beam %d %g %g %g\n", iBeam, alpha, C, S);
    
      fullMatrix<double> R(6, 6);
      R(0, 0) = R(1, 1) = R(3, 3) = R(4, 4) = C;
      R(0, 1) = R(3, 4) = S;
      R(1, 0) = R(4, 3) = -S;
      R(2, 2) = R(5, 5) = 1.0;
    
      fullMatrix<double> k(6, 6);
    
      // tensile stiffness
      k(0, 0) = k(3, 3) = TS;
      k(0, 3) = k(3, 0) = -TS;
      // bending stiffness
      k(1, 1) = k(4, 4) = 12 * BS;
      k(2, 2) = k(5, 5) = 4. * BS * b._l * b._l;
      k(1, 2) = k(2, 1) = k(1, 5) = k(5, 1) = 6 * BS * b._l;
      k(4, 2) = k(2, 4) = k(4, 5) = k(5, 4) = -6 * BS * b._l;
      k(4, 1) = k(1, 4) = -12 * BS;
      k(5, 2) = k(2, 5) = -2 * BS * b._l * b._l;
    
      fullMatrix<double> Rt(R), temp(6, 6);
      Rt.transposeInPlace();
      Rt.mult(k, temp);
      temp.mult(R, K);
    }
    
    void frameSolver2d::solve()
    {
    #if defined(HAVE_PETSC)
      linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
    #elif defined(HAVE_GMM)
      linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
      lsys->setGmres(1);
      lsys->setNoisy(1);
    #else
      linearSystemFull<double> *lsys = new linearSystemFull<double>;
    #endif
    
      if(pAssembler) delete pAssembler;
      pAssembler = new dofManager<double>(lsys);
    
      // fix dofs and create free ones
      createDofs();
    
      // force vector
      std::vector<std::pair<GVertex *, std::vector<double> > >::iterator it =
        _nodalForces.begin();
      for(; it != _nodalForces.end(); ++it) {
        MVertex *v = it->first->mesh_vertices[0];
        const std::vector<double> &F = it->second;
        Dof DOFX(v->getNum(), 0);
        Dof DOFY(v->getNum(), 1);
        pAssembler->assemble(DOFX, F[0]);
        pAssembler->assemble(DOFY, F[1]);
      }
    
      // stifness matrix
      for(std::size_t i = 0; i < _beams.size(); i++) {
        fullMatrix<double> K(6, 6);
        computeStiffnessMatrix(i, K);
        _beams[i]._stiffness = K;
        MVertex *v0 = _beams[i]._element->getVertex(0);
        MVertex *v1 = _beams[i]._element->getVertex(1);
        Dof theta0(v0->getNum(),
                   Dof::createTypeWithTwoInts(2, _beams[i]._rotationTags[0]));
        Dof theta1(v1->getNum(),
                   Dof::createTypeWithTwoInts(2, _beams[i]._rotationTags[1]));
        Dof U0(v0->getNum(), 0);
        Dof U1(v1->getNum(), 0);
        Dof V0(v0->getNum(), 1);
        Dof V1(v1->getNum(), 1);
        Dof DOFS[6] = {U0, V0, theta0, U1, V1, theta1};
        for(int j = 0; j < 6; j++) {
          for(int k = 0; k < 6; k++) {
            pAssembler->assemble(DOFS[j], DOFS[k], K(j, k));
          }
        }
      }
      lsys->systemSolve();
    
      // save the solution
      for(std::size_t i = 0; i < _beams.size(); i++) {
        MVertex *v0 = _beams[i]._element->getVertex(0);
        MVertex *v1 = _beams[i]._element->getVertex(1);
        Dof theta0(v0->getNum(),
                   Dof::createTypeWithTwoInts(2, _beams[i]._rotationTags[0]));
        Dof theta1(v1->getNum(),
                   Dof::createTypeWithTwoInts(2, _beams[i]._rotationTags[1]));
        Dof U0(v0->getNum(), 0);
        Dof U1(v1->getNum(), 0);
        Dof V0(v0->getNum(), 1);
        Dof V1(v1->getNum(), 1);
        Dof DOFS[6] = {U0, V0, theta0, U1, V1, theta1};
        for(int j = 0; j < 6; j++) {
          pAssembler->getDofValue(DOFS[j], _beams[i]._displacement[j]);
        }
      }
      delete lsys;
      delete pAssembler;
    }
    
    void frameSolver2d::exportFrameData(const char *DISPL, const char *M)
    {
    #if defined(HAVE_POST)
      {
        std::map<int, std::vector<double> > data;
        for(std::size_t i = 0; i < _beams.size(); i++) {
          std::vector<double> tmp;
          // tmp.push_back(_beams[i]._e);
          // tmp.push_back(_beams[i]._i);
          // tmp.push_back(_beams[i]._a);
          tmp.reserve(6);
          for(int j = 0; j < 6; j++) {
            tmp.push_back(_beams[i]._displacement[j]);
          }
          data[_beams[i]._element->getNum()] = tmp;
        }
        PView *pv = new PView("displacements", "Beam", _myModel, data, 0.0, 6);
        pv->getData()->writeMSH(DISPL);
        delete pv;
      }
      {
        std::map<int, std::vector<double> > data;
        for(std::size_t i = 0; i < _beams.size(); i++) {
          std::vector<double> tmp;
          fullVector<double> d(_beams[i]._displacement, 6), F(6);
          _beams[i]._stiffness.mult(d, F);
          tmp.push_back(-F(2));
          tmp.push_back(F(5));
          data[_beams[i]._element->getNum()] = tmp;
        }
        PView *pv =
          new PView("Momentum", "ElementNodeData", _myModel, data, 0.0, 1);
        pv->getData()->writeMSH(M);
        delete pv;
      }
    #endif
    }
    
    void frameSolver2d::computeRotationTags()
    {
      std::multimap<MVertex *, gmshBeam2d *> v2b;
      for(std::size_t i = 0; i < _beams.size(); i++) {
        v2b.insert(std::make_pair(_beams[i]._element->getVertex(0), &_beams[i]));
        v2b.insert(std::make_pair(_beams[i]._element->getVertex(1), &_beams[i]));
      }
    
      std::multimap<MVertex *, gmshBeam2d *>::iterator s_it;
      for(std::multimap<MVertex *, gmshBeam2d *>::iterator it = v2b.begin();
          it != v2b.end(); it = s_it) {
        MVertex *theKey = it->first;
    
        std::pair<std::multimap<MVertex *, gmshBeam2d *>::iterator,
                  std::multimap<MVertex *, gmshBeam2d *>::iterator>
          keyRange = v2b.equal_range(theKey);
        int countRotules = 0;
        for(s_it = keyRange.first; s_it != keyRange.second; ++s_it) {
          gmshBeam2d *b = s_it->second;
          if(!b->isRigid(theKey)) {
            b->setRotationTag(theKey, ++countRotules);
          }
        }
      }
    }