Skip to content
Snippets Groups Projects
Select Git revision
  • gmsh_4_7_1
  • master default protected
  • patch_releases_4_14
  • overlaps_tags_and_distributed_export
  • overlaps_tags_and_distributed_export_rebased
  • relaying
  • alphashapes
  • 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

t20.geo

Blame
  • MEdge.cpp 5.82 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 <algorithm>
    #include "MEdge.h"
    #include "Numeric.h"
    #include "GmshDefines.h"
    #include "ElementType.h"
    #include "nodalBasis.h"
    #include "BasisFactory.h"
    
    // FIXME
    // remove that when MElementCut is removed
    bool MEdge::isInside(MVertex *v) const
    {
      double tol = MVertexLessThanLexicographic::getTolerance();
      MVertex *v0 = _v[0];
      MVertex *v1 = _v[1];
      MVertexLessThanLexicographic lt;
      if(lt(v0, v1)) {
        v0 = _v[1];
        v1 = _v[0];
      }
      double x = v->x(), y = v->y(), z = v->z();
      double x0 = v0->x(), y0 = v0->y(), z0 = v0->z();
      double x1 = v1->x(), y1 = v1->y(), z1 = v1->z();
      if(fabs(x - x0) < tol && fabs(y - y0) < tol && fabs(z - z0) < tol)
        return true;
      if(fabs(x - x1) < tol && fabs(y - y1) < tol && fabs(z - z1) < tol)
        return true;
      if(x < x0 - tol || x > x1 + tol || y < std::min(y0, y1) - tol ||
         y > std::max(y0, y1) + tol || z < std::min(z0, z1) - tol ||
         z > std::max(z0, z1) + tol)
        return false;
      if(fabs(x1 - x0) > tol) {
        double tx = (x - x0) / (x1 - x0);
        if(fabs(y1 - y0) > tol) {
          double ty = (y - y0) / (y1 - y0);
          if(fabs(z1 - z0) > tol) {
            double tz = (z - z0) / (z1 - z0);
            if(fabs(tx - ty) > tol || fabs(tx - tz) > tol) return false;
          }
          else {
            if(fabs(tx - ty) > tol) return false;
          }
        }
        else {
          if(fabs(z1 - z0) > tol) {
            double tz = (z - z0) / (z1 - z0);
            if(fabs(tx - tz) > tol) return false;
          }
        }
      }
      else {
        if(fabs(y1 - y0) > tol) {
          double ty = (y - y0) / (y1 - y0);
          if(fabs(z1 - z0) > tol) {
            double tz = (z - z0) / (z1 - z0);
            if(fabs(ty - tz) > tol) return false;
          }
        }
      }
      return true;
    }
    
    bool SortEdgeConsecutive(const std::vector<MEdge> &e,
                             std::vector<std::vector<MVertex *> > &vs)
    {
      if(e.empty()) return true;
      std::map<MVertex *, std::pair<MVertex *, MVertex *> > c;
    
      for(size_t i = 0; i < e.size(); i++) {
        MVertex *v0 = e[i].getVertex(0);
        MVertex *v1 = e[i].getVertex(1);
        std::map<MVertex *, std::pair<MVertex *, MVertex *> >::iterator it0 =
          c.find(v0);
        std::map<MVertex *, std::pair<MVertex *, MVertex *> >::iterator it1 =
          c.find(v1);
        if(it0 == c.end())
          c[v0] = std::make_pair(v1, (MVertex *)NULL);
        else {
          if(it0->second.second == NULL) { it0->second.second = v1; }
          else {
            Msg::Debug("A list of edges has points that are adjacent to 3 edges");
            return false;
          }
        }
        if(it1 == c.end())
          c[v1] = std::make_pair(v0, (MVertex *)NULL);
        else {
          if(it1->second.second == NULL) { it1->second.second = v0; }
          else {
            Msg::Debug("Wrong topology for a list of edges");
            Msg::Debug("Node %d is adjacent to more than 2 nodes %d %d",
                       v1->getNum(),it1->second.first->getNum(),
                       it1->second.second->getNum());
            return false;
          }
        }
      }
    
      while(!c.empty()) {
        std::vector<MVertex *> v;
        MVertex *start = NULL;
        {
          std::map<MVertex *, std::pair<MVertex *, MVertex *> >::iterator it =
            c.begin();
          start = it->first;
          for(; it != c.end(); ++it) {
            if(it->second.second == NULL) {
              start = it->first;
              break;
            }
          }
        }
    
        std::map<MVertex *, std::pair<MVertex *, MVertex *> >::iterator its =
          c.find(start);
    
        MVertex *prev =
          (its->second.second == start) ? its->second.first : its->second.second;
        MVertex *current = start;
    
        do {
          if(c.size() == 0) {
            Msg::Warning("Wrong topology in a wire");
            return false;
          }
          v.push_back(current);
          std::map<MVertex *, std::pair<MVertex *, MVertex *> >::iterator it =
            c.find(current);
          if(it == c.end() || it->first == NULL) {
            Msg::Error("Impossible to find %d", current->getNum());
            return false;
          }
          MVertex *v1 = it->second.first;
          MVertex *v2 = it->second.second;
          c.erase(it);
          MVertex *temp = current;
          if(v1 == prev)
            current = v2;
          else if(v2 == prev)
            current = v1;
          else {
            break;
          }
          prev = temp;
          if(current == start) { v.push_back(current); }
        } while(current != start && current != NULL);
        vs.push_back(v);
      }
      return true;
    }
    
    MEdgeN::MEdgeN(const std::vector<MVertex *> &v)
    {
      _v.resize(v.size());
      for(std::size_t i = 0; i < v.size(); i++) _v[i] = v[i];
    }
    
    MEdge MEdgeN::getEdge() const { return MEdge(_v[0], _v[1]); }
    
    SPoint3 MEdgeN::pnt(double u) const
    {
      int tagLine = ElementType::getType(TYPE_LIN, getPolynomialOrder());
      const nodalBasis *fs = BasisFactory::getNodalBasis(tagLine);
    
      double f[100];
      fs->f(u, 0, 0, f);
    
      double x = 0, y = 0, z = 0;
      for(int i = 0; i < fs->getNumShapeFunctions(); i++) {
        x += f[i] * _v[i]->x();
        y += f[i] * _v[i]->y();
        z += f[i] * _v[i]->z();
      }
      return SPoint3(x, y, z);
    }
    
    SVector3 MEdgeN::tangent(double u) const
    {
      int tagLine = ElementType::getType(TYPE_LIN, getPolynomialOrder());
      const nodalBasis *fs = BasisFactory::getNodalBasis(tagLine);
    
      double sf[100][3];
      fs->df(u, 0, 0, sf);
    
      double dx = 0, dy = 0, dz = 0;
      for(int i = 0; i < fs->getNumShapeFunctions(); i++) {
        dx += sf[i][0] * _v[i]->x();
        dy += sf[i][0] * _v[i]->y();
        dz += sf[i][0] * _v[i]->z();
      }
      return SVector3(dx, dy, dz).unit();
    }
    
    double MEdgeN::interpolate(const double val[], double u, int stride) const
    {
      int tagLine = ElementType::getType(TYPE_LIN, getPolynomialOrder());
      const nodalBasis *fs = BasisFactory::getNodalBasis(tagLine);
    
      double f[100];
      fs->f(u, 0, 0, f);
    
      double sum = 0;
      for(int i = 0, k = 0; i < fs->getNumShapeFunctions(); i++, k += stride) {
        sum += f[i] * val[k];
      }
      return sum;
    }