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

t6.jl

Blame
  • discreteEdge.cpp 13.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 <vector>
    #include "GmshConfig.h"
    #include "GmshMessage.h"
    #include "discreteEdge.h"
    #include "MLine.h"
    #include "MPoint.h"
    #include "GModelIO_GEO.h"
    #include "Geo.h"
    
    #if defined(HAVE_MESH)
    #include "meshGEdge.h"
    #include "Context.h"
    #endif
    
    discreteEdge::discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1)
      : GEdge(model, num, _v0, _v1)
    {
      Curve *c = CreateCurve(num, MSH_SEGM_DISCRETE, 0, 0, 0, -1, -1, 0., 1.);
      Tree_Add(model->getGEOInternals()->Curves, &c);
      CreateReversedCurve(c);
      _split[0] = _split[1] = NULL;
    }
    
    discreteEdge::discreteEdge(GModel *model, int num) : GEdge(model, num)
    {
      Curve *c = CreateCurve(num, MSH_SEGM_DISCRETE, 0, 0, 0, -1, -1, 0., 1.);
      Tree_Add(model->getGEOInternals()->Curves, &c);
      CreateReversedCurve(c);
      _split[0] = _split[1] = NULL;
    }
    
    discreteEdge::~discreteEdge()
    {
      _split[0] = _split[1] = NULL;
    }
    
    void discreteEdge::orderMLines()
    {
      size_t ss = lines.size();
      if(!ss) return;
    
      std::vector<MEdge> ed;
      std::vector<std::vector<MVertex *> > vs;
      for(std::size_t i = 0; i < lines.size(); i++) {
        ed.push_back(MEdge(lines[i]->getVertex(0), lines[i]->getVertex(1)));
        delete lines[i];
      }
      lines.clear();
    
      if(!SortEdgeConsecutive(ed, vs))
        Msg::Warning("Discrete curve elements cannot be ordered");
    
      if(vs.size() != 1)
        Msg::Warning("Discrete curve %d is mutiply connected", tag());
    
      std::size_t START = 0;
      for(; START < vs[0].size(); START++)
        if(vs[0][START]->onWhat()->dim() == 0) break;
    
      if(START == vs[0].size())
        Msg::Warning("Discrete curve %d topology is wrong", tag());
    
      std::size_t i = START;
      while(lines.size() != ss) {
        if(vs[0][i % vs[0].size()] != vs[0][(i + 1) % vs[0].size()])
          lines.push_back(
            new MLine(vs[0][i % vs[0].size()], vs[0][(i + 1) % vs[0].size()]));
        i++;
      }
    
      mesh_vertices.clear();
      for(std::size_t i = 0; i < lines.size() - 1; ++i) {
        MVertex *v11 = lines[i]->getVertex(1);
        mesh_vertices.push_back(v11);
      }
      if(lines.empty()) {
        Msg::Error("No line elements in discrete curve %d", tag());
        return;
      }
      GVertex *g0 = dynamic_cast<GVertex *>(lines[0]->getVertex(0)->onWhat());
      if(g0) setBeginVertex(g0);
      GVertex *g1 =
        dynamic_cast<GVertex *>(lines[lines.size() - 1]->getVertex(1)->onWhat());
      if(g1) setEndVertex(g1);
      if(!g0 || !g1)
        Msg::Error("Discrete curve %d has non consecutive line elements", tag());
    }
    
    bool discreteEdge::getLocalParameter(const double &t, int &iLine,
                                         double &tLoc) const
    {
      for(iLine = 0; iLine < (int)_discretization.size() - 1; iLine++) {
        double tmin = _pars[iLine];
        double tmax = _pars[iLine + 1];
        if(t >= tmin && t <= tmax) {
          tLoc = (t - tmin) / (tmax - tmin);
          return true;
        }
      }
      return false;
    }
    
    GPoint discreteEdge::point(double par) const
    {
      double tLoc;
      int iEdge;
    
      if(!getLocalParameter(par, iEdge, tLoc)) return GPoint();
    
      SPoint3 vB = _discretization[iEdge];
      SPoint3 vE = _discretization[iEdge + 1];
    
      // linear Lagrange mesh
      SPoint3 v = vB + (vE - vB) * tLoc;
      return GPoint(v.x(), v.y(), v.z(), this, par);
    }
    
    SVector3 discreteEdge::firstDer(double par) const
    {
      double tLoc;
      int iEdge;
    
      if(!getLocalParameter(par, iEdge, tLoc)) return SVector3();
    
      SPoint3 vB = _discretization[iEdge];
      SPoint3 vE = _discretization[iEdge + 1];
    
      double dx, dy, dz;
      dx = (vE.x() - vB.x()); // / dt;
      dy = (vE.y() - vB.y()); // / dt;
      dz = (vE.z() - vB.z()); // / dt;
    
      SVector3 der(dx, dy, dz);
      return der;
    }
    
    SPoint2 discreteEdge::reparamOnFace(const GFace *face, double epar,
                                        int dir) const
    {
      GPoint p = point(epar);
      double guess[2];
      GPoint ps = face->closestPoint(SPoint3(p.x(), p.y(), p.z()), guess);
      return SPoint2(ps.u(), ps.v());
    }
    
    double discreteEdge::curvature(double par) const
    {
      return 1.e-12;
      double tLoc;
      int iEdge;
      if(_discretization.size() <= 3)
        return 0.0; // no clue on how to compute curvature with so few data...
    
      if(!getLocalParameter(par, iEdge, tLoc)) return 0.0;
    
      // Take 3 points x y z : radius of curvature is equal to  |x-y| |x-z| |z-y| /
      // area (x,y,z)
    
      int iEdgePlus = iEdge + 1, iEdgeMinus = iEdge - 1;
    
      if(iEdge == 0) {
        iEdge++;
        iEdgePlus++;
        iEdgeMinus++;
        if(periodic(0)) { iEdgeMinus = _discretization.size() - 2; }
      }
      else if(iEdge == (int)(_discretization.size() - 2)) {
        iEdge--;
        iEdgePlus--;
        iEdgeMinus--;
        if(periodic(0)) { iEdgePlus = 0; }
      }
    
      SPoint3 a =
        (_discretization[iEdgeMinus] + _discretization[iEdgeMinus + 1]) * .5;
      SPoint3 b = (_discretization[iEdge] + _discretization[iEdge + 1]) * .5;
      SPoint3 c =
        (_discretization[iEdgePlus] + _discretization[iEdgePlus + 1]) * .5;
    
      // SPoint3 a((discrete_lines[iEdgeMinus]->getVertex(0)->x() +
      //            discrete_lines[iEdgeMinus]->getVertex(1)->x()) *
      //             .5,
      //           (discrete_lines[iEdgeMinus]->getVertex(0)->y() +
      //            discrete_lines[iEdgeMinus]->getVertex(1)->y()) *
      //             .5,
      //           (discrete_lines[iEdgeMinus]->getVertex(0)->z() +
      //            discrete_lines[iEdgeMinus]->getVertex(1)->z()) *
      //             .5);
      // SPoint3 b((discrete_lines[iEdge]->getVertex(0)->x() +
      //            discrete_lines[iEdge]->getVertex(1)->x()) *
      //             .5,
      //           (discrete_lines[iEdge]->getVertex(0)->y() +
      //            discrete_lines[iEdge]->getVertex(1)->y()) *
      //             .5,
      //           (discrete_lines[iEdge]->getVertex(0)->z() +
      //            discrete_lines[iEdge]->getVertex(1)->z()) *
      //             .5);
      // SPoint3 c((discrete_lines[iEdgePlus]->getVertex(0)->x() +
      //            discrete_lines[iEdgePlus]->getVertex(1)->x()) *
      //             .5,
      //           (discrete_lines[iEdgePlus]->getVertex(0)->y() +
      //            discrete_lines[iEdgePlus]->getVertex(1)->y()) *
      //             .5,
      //           (discrete_lines[iEdgePlus]->getVertex(0)->z() +
      //            discrete_lines[iEdgePlus]->getVertex(1)->z()) *
      //             .5);
    
      double A = b.distance(c);
      double B = c.distance(a);
      double C = a.distance(b);
    
      // radius of the circumcircle ...
      // Heron's formula for the area of a triangle
      double R =
        A * B * C / sqrt((A + B + C) * (-A + B + C) * (A - B + C) * (A + B - C));
    
      // printf("R(%d,%g) = %g\n",tag(),par,R);
    
      return R = 0.0 ? 1.E22 : 1. / R;
    }
    
    Range<double> discreteEdge::parBounds(int i) const
    {
      return Range<double>(0, (double)(_discretization.size() - 1));
    }
    
    void discreteEdge::createGeometry()
    {
      std::vector<MVertex *> discrete_vertices;
      std::vector<MLine *> discrete_lines;
    
      if(lines.empty()) return;
    
      if(!_discretization.empty()) return;
    
      orderMLines();
    
      bool reverse = false;
      if(getEndVertex())
        reverse = (lines[0]->getVertex(0) == getEndVertex()->mesh_vertices[0]);
    
      std::map<MVertex *, MVertex *> old2new;
    
      for(std::size_t i = 0; i < lines.size(); i++) {
        for(std::size_t j = 0; j < 2; j++) {
          MVertex *v = lines[i]->getVertex(j);
          if(old2new.find(v) == old2new.end()) {
            MVertex *vnew = new MVertex(v->x(), v->y(), v->z(), this);
            old2new[v] = vnew;
          }
        }
      }
    
      std::vector<MLine *> _temp;
      discrete_lines.resize(lines.size());
      for(std::size_t i = 0; i < lines.size(); i++) {
        MVertex *v0 = lines[i]->getVertex(0);
        MVertex *v1 = lines[i]->getVertex(1);
        MVertex *v00 = old2new[v0];
        MVertex *v01 = old2new[v1];
        if(reverse)
          discrete_lines[lines.size() - i - 1] = new MLine(v01, v00);
        else
          discrete_lines[i] = new MLine(v00, v01);
      }
    
      // compute parameters and recompute the vertices
      _discretization.push_back(SPoint3(discrete_lines[0]->getVertex(0)->x(),
                                        discrete_lines[0]->getVertex(0)->y(),
                                        discrete_lines[0]->getVertex(0)->z()));
    
      for(std::size_t i = 0; i < discrete_lines.size(); i++) {
        _discretization.push_back(SPoint3(discrete_lines[i]->getVertex(1)->x(),
                                          discrete_lines[i]->getVertex(1)->y(),
                                          discrete_lines[i]->getVertex(1)->z()));
      }
    
      _pars.push_back(0.0);
      for(std::size_t i = 1; i < discrete_lines.size(); i++) {
        _pars.push_back((double)i);
        MVertex *newv = discrete_lines[i]->getVertex(0);
        discrete_vertices.push_back(newv);
      }
      _pars.push_back((double)discrete_lines.size());
      //  mesh_vertices.clear();
      //  lines.clear();
      for(std::size_t i = 0; i < discrete_lines.size(); i++)
        delete discrete_lines[i];
      for(std::size_t i = 0; i < discrete_vertices.size(); i++)
        delete discrete_vertices[i];
      discrete_lines.clear();
      discrete_vertices.clear();
    }
    
    void discreteEdge::mesh(bool verbose)
    {
    #if defined(HAVE_MESH)
      if(_discretization.empty()) return;
      meshGEdge mesher;
      mesher(this);
    #endif
    }
    
    bool discreteEdge::split(MVertex *v, GVertex *gv, int &TAG)
    {
      GVertex *gv0 = getBeginVertex();
      GVertex *gv1 = getEndVertex();
    
      if(v != gv->mesh_vertices[0]) {
        Msg::Error("Wrong vertex for splitting discrete curve");
        return false;
      }
      discreteEdge *de_new[2];
    
      de_new[0] = new discreteEdge(model(), ++TAG, gv0, gv);
      de_new[1] = new discreteEdge(model(), ++TAG, gv, gv1);
    
      setSplit(de_new[0], de_new[1], gv);
    
      int current = 0;
      de_new[current]->lines.push_back(lines[0]);
      for(size_t i = 1; i < lines.size(); i++) {
        if(lines[i]->getVertex(0) == v) { current++; }
        else {
          de_new[current]->mesh_vertices.push_back(lines[i]->getVertex(0));
          lines[i]->getVertex(0)->setEntity(de_new[current]);
        }
        de_new[current]->lines.push_back(lines[i]);
      }
      lines.clear();
      mesh_vertices.clear();
      Msg::Info("Edge %d split in %d %d (%d and %d lines)\n", tag(),
                de_new[0]->tag(), de_new[1]->tag(), de_new[0]->lines.size(),
                de_new[1]->lines.size());
    
      std::vector<GFace *> f = faces();
      for(size_t i = 0; i < f.size(); i++) {
        std::vector<GEdge *> new_eds, old_eds;
        old_eds = f[i]->edges();
        for(size_t j = 0; j < old_eds.size(); j++) {
          if(old_eds[j] == this) {
            new_eds.push_back(de_new[0]);
            new_eds.push_back(de_new[1]);
            new_eds[0]->addFace(f[i]);
            new_eds[1]->addFace(f[i]);
          }
          else
            new_eds.push_back(old_eds[j]);
        }
        f[i]->set(new_eds);
      }
      model()->add(de_new[0]);
      model()->add(de_new[1]);
      model()->remove(this);
      return true;
    }
    
    void discreteEdge::unsplit(void)
    {
      if(_split[2] != NULL) return;
      if(_split[0] == NULL) return;
      std::vector<GEdge *> l_edges;
      std::vector<GVertex *> l_vertices;
      getSplit(l_edges, l_vertices);
      // remove all internal model vertices
      for(size_t k = 0; k < l_vertices.size(); k++) {
        l_vertices[k]->mesh_vertices.clear();
        model()->remove(l_vertices[k]);
      }
      // remove all internal model edges
      for(size_t k = 0; k < l_edges.size(); k++) {
        for(size_t l = 0; l < l_edges[k]->lines.size(); l++) {
          lines.push_back(l_edges[k]->lines[l]);
          if(l_edges[k]->lines[l]->getVertex(0)->onWhat() != getBeginVertex() &&
             l_edges[k]->lines[l]->getVertex(0)->onWhat() != getEndVertex() &&
             l_edges[k]->lines[l]->getVertex(0)->onWhat() != this)
            l_edges[k]->lines[l]->getVertex(0)->setEntity(this);
          if(l_edges[k]->lines[l]->getVertex(1)->onWhat() != getBeginVertex() &&
             l_edges[k]->lines[l]->getVertex(1)->onWhat() != getEndVertex() &&
             l_edges[k]->lines[l]->getVertex(1)->onWhat() != this)
            l_edges[k]->lines[l]->getVertex(1)->setEntity(this);
          if(l_edges[k]->lines[l]->getVertex(1)->onWhat() == this &&
             std::find(mesh_vertices.begin(), mesh_vertices.end(),
                       l_edges[k]->lines[l]->getVertex(0)) == mesh_vertices.end())
            mesh_vertices.push_back(l_edges[k]->lines[l]->getVertex(1));
          if(l_edges[k]->lines[l]->getVertex(1)->onWhat() == this &&
             std::find(mesh_vertices.begin(), mesh_vertices.end(),
                       l_edges[k]->lines[l]->getVertex(1)) == mesh_vertices.end())
            mesh_vertices.push_back(l_edges[k]->lines[l]->getVertex(1));
        }
        l_edges[k]->lines.clear();
        l_edges[k]->mesh_vertices.clear();
        model()->remove(l_edges[k]);
      }
      // recompute the right edge 2 face topology
    
      std::vector<GFace *> f = faces();
      for(size_t i = 0; i < f.size(); i++) {
        std::vector<GEdge *> new_eds, old_eds;
        old_eds = f[i]->edges();
        new_eds.push_back(this);
        for(size_t j = 0; j < old_eds.size(); j++) {
          if(std::find(l_edges.begin(), l_edges.end(), old_eds[j]) ==
             l_edges.end()) {
            new_eds.push_back(old_eds[j]);
          }
        }
        f[i]->set(new_eds);
      }
    }
    
    void discreteEdge::writeParametrization(FILE *fp, bool binary)
    {
      fprintf(fp, "%lu\n", _discretization.size());
      for(size_t i = 0; i < _discretization.size(); i++) {
        fprintf(fp, "%22.15E %22.15E %22.15E %22.15E\n", _discretization[i].x(), _discretization[i].y(),
                _discretization[i].z(), _pars[i]);
      }
    }
    
    void discreteEdge::readParametrization(FILE *fp, bool binary)
    {
      int N;
      if(fscanf(fp, "%d", &N) != 1)
        return;
    
      _pars.resize(N);
      _discretization.resize(N);
      for(int i = 0; i < N; i++) {
        double x, y, z, t;
        if(fscanf(fp, "%lf %lf %lf %lf\n", &x, &y, &z, &t) != 4) {
          return;
        }
        _pars[i] = t;
        _discretization[i] = SPoint3(x, y, z);
      }
    }