Skip to content
Snippets Groups Projects
Select Git revision
  • 427de34d56536a6c7588383b862743eb16df200f
  • master default protected
  • pluginMeshQuality
  • fixBugsAmaury
  • 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
  • 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

Navigator.cpp

Blame
  • Crack.cpp 9.13 KiB
    // Gmsh - Copyright (C) 1997-2018 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file for license information. Please report all
    // bugs and problems to the public mailing list <gmsh@onelab.info>.
    
    #include "Crack.h"
    #include "GModel.h"
    #include "discreteEdge.h"
    #include "discreteFace.h"
    #include "MElement.h"
    #include "MLine.h"
    #include "MTriangle.h"
    #include "MQuadrangle.h"
    #include "MEdge.h"
    #include "Context.h"
    
    StringXNumber CrackOptions_Number[] = {
      {GMSH_FULLRC, "Dimension", NULL, 1.},
      {GMSH_FULLRC, "PhysicalGroup", NULL, 1.},
      {GMSH_FULLRC, "OpenBoundaryPhysicalGroup", NULL, 0.},
      {GMSH_FULLRC, "NormalX", NULL, 0.},
      {GMSH_FULLRC, "NormalY", NULL, 0.},
      {GMSH_FULLRC, "NormalZ", NULL, 1.}};
    
    extern "C" {
    GMSH_Plugin *GMSH_RegisterCrackPlugin() { return new GMSH_CrackPlugin(); }
    }
    
    std::string GMSH_CrackPlugin::getHelp() const
    {
      return "Plugin(Crack) creates a crack around the physical "
             "group `PhysicalGroup' of dimension `Dimension' (1 or 2), "
             "embedded in a mesh of dimension `Dimension' + 1. "
             "The plugin duplicates the vertices and the elements on "
             "the crack and stores them in a new discrete curve "
             "(`Dimension' = 1) or surface (`Dimension' = 2). The "
             "elements touching the crack on the ``negative'' side "
             "are modified to use the newly generated vertices."
             "If `OpenBoundaryPhysicalGroup' is given (> 0), its "
             "vertices are duplicated and the crack will be left "
             "open on that (part of the) boundary. Otherwise, the "
             "lips of the crack are sealed, i.e., its vertices are "
             "not duplicated. For 1D cracks, `NormalX', `NormalY' and "
             "`NormalZ' provide the reference normal of the surface "
             "in which the crack is supposed to be embedded.";
    }
    
    int GMSH_CrackPlugin::getNbOptions() const
    {
      return sizeof(CrackOptions_Number) / sizeof(StringXNumber);
    }
    
    StringXNumber *GMSH_CrackPlugin::getOption(int iopt)
    {
      return &CrackOptions_Number[iopt];
    }
    
    class EdgeData {
    public:
      EdgeData(MEdge e) : edge(e) {}
      MEdge edge;
      std::vector<MVertex *> data;
    };
    
    struct Less_EdgeData : public std::binary_function<EdgeData, EdgeData, bool> {
      bool operator()(const EdgeData &e1, const EdgeData &e2) const
      {
        if(e1.edge.getMinVertex() < e2.edge.getMinVertex()) return true;
        if(e1.edge.getMinVertex() > e2.edge.getMinVertex()) return false;
        if(e1.edge.getMaxVertex() < e2.edge.getMaxVertex()) return true;
        return false;
      }
    };
    
    PView *GMSH_CrackPlugin::execute(PView *view)
    {
      int dim = (int)CrackOptions_Number[0].def;
      int physical = (int)CrackOptions_Number[1].def;
      int open = (int)CrackOptions_Number[2].def;
      SVector3 normal1d(CrackOptions_Number[3].def, CrackOptions_Number[4].def,
                        CrackOptions_Number[5].def);
    
      if(dim != 1 && dim != 2) {
        Msg::Error("Crack dimension should be 1 or 2");
        return view;
      }
    
      GModel *m = GModel::current();
      std::map<int, std::vector<GEntity *> > groups[4];
      m->getPhysicalGroups(groups);
      std::vector<GEntity *> entities = groups[dim][physical];
    
      if(entities.empty()) {
        Msg::Error("Physical group %d (dimension %d) is empty", physical, dim);
        return view;
      }
    
      std::vector<GEntity *> openEntities;
      if(open > 0) {
        openEntities = groups[dim - 1][open];
        if(openEntities.empty()) {
          Msg::Error("Open boundary physical group %d (dimension %d) is empty",
                     open, dim - 1);
          return view;
        }
      }
    
      std::set<GEntity *> crackEntities;
      crackEntities.insert(entities.begin(), entities.end());
      crackEntities.insert(openEntities.begin(), openEntities.end());
    
      // get crack elements
      std::vector<MElement *> crackElements;
      for(unsigned int i = 0; i < entities.size(); i++)
        for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++)
          crackElements.push_back(entities[i]->getMeshElement(j));
    
      // get internal crack vertices and boundary vertices
      std::set<MVertex *> crackVertices, bndVertices;
      if(dim == 1) {
        for(unsigned int i = 0; i < crackElements.size(); i++) {
          for(std::size_t j = 0; j < crackElements[i]->getNumVertices(); j++) {
            MVertex *v = crackElements[i]->getVertex(j);
            crackVertices.insert(v);
          }
          for(int j = 0; j < crackElements[i]->getNumPrimaryVertices(); j++) {
            MVertex *v = crackElements[i]->getVertex(j);
            if(bndVertices.find(v) == bndVertices.end())
              bndVertices.insert(v);
            else
              bndVertices.erase(v);
          }
        }
      }
      else {
        std::set<EdgeData, Less_EdgeData> bnd;
        for(unsigned int i = 0; i < crackElements.size(); i++) {
          for(std::size_t j = 0; j < crackElements[i]->getNumVertices(); j++) {
            MVertex *v = crackElements[i]->getVertex(j);
            crackVertices.insert(v);
          }
          for(int j = 0; j < crackElements[i]->getNumEdges(); j++) {
            EdgeData ed(crackElements[i]->getEdge(j));
            if(bnd.find(ed) == bnd.end()) {
              crackElements[i]->getEdgeVertices(j, ed.data);
              bnd.insert(ed);
            }
            else
              bnd.erase(ed);
          }
        }
        for(std::set<EdgeData, Less_EdgeData>::iterator it = bnd.begin();
            it != bnd.end(); it++)
          bndVertices.insert(it->data.begin(), it->data.end());
      }
    
      // get (forced) open boundary vertices and remove them from boundary vertices
      for(unsigned int i = 0; i < openEntities.size(); i++) {
        for(unsigned int j = 0; j < openEntities[i]->getNumMeshElements(); j++) {
          MElement *e = openEntities[i]->getMeshElement(j);
          for(std::size_t k = 0; k < e->getNumVertices(); k++) {
            MVertex *v = e->getVertex(k);
            bndVertices.erase(v);
          }
        }
      }
      for(std::set<MVertex *>::iterator it = bndVertices.begin();
          it != bndVertices.end(); it++)
        crackVertices.erase(*it);
    
      // compute elements on one side of the crack
      std::set<MElement *> oneside;
      std::vector<GEntity *> allentities;
      m->getEntities(allentities);
      for(unsigned int ent = 0; ent < allentities.size(); ent++) {
        if(crackEntities.find(allentities[ent]) != crackEntities.end()) continue;
        for(unsigned int i = 0; i < allentities[ent]->getNumMeshElements(); i++) {
          MElement *e = allentities[ent]->getMeshElement(i);
          for(std::size_t j = 0; j < e->getNumVertices(); j++) {
            if(crackVertices.find(e->getVertex(j)) != crackVertices.end()) {
              // element touches the crack: find the closest crack element
              SPoint3 b = e->barycenter();
              double d = 1e200;
              MElement *ce = 0;
              for(unsigned int k = 0; k < crackElements.size(); k++) {
                double d2 = b.distance(crackElements[k]->barycenter());
                if(d2 < d) {
                  d = d2;
                  ce = crackElements[k];
                }
              }
              SVector3 dv = SVector3(e->barycenter(), ce->barycenter());
              SVector3 n;
              if(dim == 1)
                n = crossprod(normal1d, ce->getEdge(0).tangent());
              else
                n = ce->getFace(0).normal();
              if(dot(n, dv) < 0) {
                oneside.insert(e);
              }
            }
          }
        }
      }
    
      /*
      FILE *fp = fopen("debug.pos", "w");
      if(fp){
        fprintf(fp, "View \"Ele < 0\" {\n");
        for(std::set<MElement*>::iterator it = oneside.begin(); it != oneside.end();
      it++)
          (*it)->writePOS(fp, false, true, false, false, false, false);
        fprintf(fp, "};\n");
        fclose(fp);
      }
      */
    
      // create new crack entity
      GEdge *crackEdge = 0;
      GFace *crackFace = 0;
      if(dim == 1) {
        crackEdge = new discreteEdge(m, m->getMaxElementaryNumber(1) + 1, 0, 0);
        m->add(crackEdge);
      }
      else {
        crackFace = new discreteFace(m, m->getMaxElementaryNumber(2) + 1);
        m->add(crackFace);
      }
      GEntity *crackEntity =
        crackEdge ? (GEntity *)crackEdge : (GEntity *)crackFace;
      crackEntity->physicals.push_back(physical);
    
      // duplicate internal crack vertices
      std::map<MVertex *, MVertex *> vxv;
      for(std::set<MVertex *>::iterator it = crackVertices.begin();
          it != crackVertices.end(); it++) {
        MVertex *v = *it;
        MVertex *newv = new MVertex(v->x(), v->y(), v->z(), crackEntity);
        crackEntity->mesh_vertices.push_back(newv);
        vxv[v] = newv;
      }
    
      // duplicate crack elements
      for(unsigned int i = 0; i < crackElements.size(); i++) {
        MElement *e = crackElements[i];
        std::vector<MVertex *> verts;
        e->getVertices(verts);
        for(unsigned int j = 0; j < verts.size(); j++) {
          if(vxv.count(verts[j])) verts[j] = vxv[verts[j]];
        }
        MElementFactory f;
        MElement *newe = f.create(e->getTypeForMSH(), verts, 0, e->getPartition());
        if(crackEdge && newe->getType() == TYPE_LIN)
          crackEdge->lines.push_back((MLine *)newe);
        else if(crackFace && newe->getType() == TYPE_TRI)
          crackFace->triangles.push_back((MTriangle *)newe);
        else if(crackFace && newe->getType() == TYPE_QUA)
          crackFace->quadrangles.push_back((MQuadrangle *)newe);
      }
    
      // replace vertices in elements on one side of the crack
      for(std::set<MElement *>::iterator it = oneside.begin(); it != oneside.end();
          it++) {
        MElement *e = *it;
        for(std::size_t i = 0; i < e->getNumVertices(); i++) {
          if(vxv.count(e->getVertex(i))) e->setVertex(i, vxv[e->getVertex(i)]);
        }
      }
    
      CTX::instance()->mesh.changed = ENT_ALL;
    
      return view;
    }