Skip to content
Snippets Groups Projects
Select Git revision
  • a5880d08663440ff07a3e78b1b723d13843957d2
  • master default protected
  • alphashapes
  • quadMeshingTools
  • cygwin_conv_path
  • macos_arm64
  • add-transfiniteautomatic-to-geo
  • patch_releases_4_10
  • HierarchicalHDiv
  • isuruf-master-patch-63355
  • hyperbolic
  • hexdom
  • hxt_update
  • jf
  • 1618-pythonocc-and-gmsh-api-integration
  • octreeSizeField
  • hexbl
  • alignIrregularVertices
  • getEdges
  • patch_releases_4_8
  • isuruf-master-patch-51992
  • 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
  • gmsh_4_8_4
  • gmsh_4_8_3
  • gmsh_4_8_2
  • gmsh_4_8_1
  • gmsh_4_8_0
  • gmsh_4_7_1
  • gmsh_4_7_0
41 results

linearSystemCSR.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    mainCartesian.cpp 11.56 KiB
    //                          lx   ly   lz    rmax  levels  refcs
    //
    // plaqueEp.stp             0.2  0.2  0.2   0.3   3
    // plaqueEpRotated.stp      0.3  0.3  0.3   0.3   3
    // jonction_collee2.stp     6    6    6     10    3
    // panneau_raidi_simple.stp 3    3    50    5     2       0
    // plaque_trouee.stp        1    1    1     2     3
    
    #include "Gmsh.h"
    #include "GModel.h"
    #include "Numeric.h"
    #include "GmshMessage.h"
    #include "cartesian.h"
    
    static void insertActiveCells(double x, double y, double z, double rmax,
                                  cartesianBox<double> &box)
    {
      int id1 = box.getCellContainingPoint(x - rmax, y - rmax, z - rmax);
      int id2 = box.getCellContainingPoint(x + rmax, y + rmax, z + rmax);
      int i1, j1 ,k1;
      box.getCellIJK(id1, i1, j1, k1);
      int i2, j2, k2;
      box.getCellIJK(id2, i2, j2, k2);
      for (int i = i1; i <= i2; i++)
        for (int j = j1; j <= j2; j++)
          for (int k = k1; k <= k2; k++)
            box.insertActiveCell(box.getCellIndex(i, j, k));
    }
    
    static void computeLevelset(GModel *gm, cartesianBox<double> &box)
    {
      // tolerance for desambiguation
      const double tol = box.getLC() * 1.e-12;
      std::vector<SPoint3> nodes;
      std::vector<int> indices;
      for (cartesianBox<double>::valIter it = box.nodalValuesBegin();
           it != box.nodalValuesEnd(); ++it){
        nodes.push_back(box.getNodeCoordinates(it->first));
        indices.push_back(it->first);
      }
      Msg::Info("  %d nodes in the grid at level %d", (int)nodes.size(), box.getLevel());
      std::vector<double> dist, localdist;
      std::vector<SPoint3> dummy;
      for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++){
        for (int i = 0; i < (*fit)->stl_triangles.size(); i += 3){
          int i1 = (*fit)->stl_triangles[i];
          int i2 = (*fit)->stl_triangles[i + 1];
          int i3 = (*fit)->stl_triangles[i + 2];
          GPoint p1 = (*fit)->point((*fit)->stl_vertices[i1]);
          GPoint p2 = (*fit)->point((*fit)->stl_vertices[i2]);
          GPoint p3 = (*fit)->point((*fit)->stl_vertices[i3]);
          SPoint2 p = ((*fit)->stl_vertices[i1] + (*fit)->stl_vertices[i2] +
                       (*fit)->stl_vertices[i3]) * 0.33333333;
          SVector3 N = (*fit)->normal(p);
          SPoint3 P1(p1.x(), p1.y(), p1.z());
          SPoint3 P2(p2.x(), p2.y(), p2.z());
          SPoint3 P3(p3.x(), p3.y(), p3.z());
          SVector3 NN(crossprod(P2 - P1, P3 - P1));
          if (dot(NN, N) > 0)
    	signedDistancesPointsTriangle(localdist, dummy, nodes, P1, P2, P3);
          else
    	signedDistancesPointsTriangle(localdist, dummy, nodes, P2, P1, P3);
          if(dist.empty())
            dist = localdist;
          else{
            for (unsigned int j = 0; j < localdist.size(); j++){
              // FIXME: if there is an ambiguity assume we are inside (to
              // avoid holes in the structure). This is definitely just a
              // hack, as it could create pockets of matter outside the
              // structure...
              if(dist[j] * localdist[j] < 0 &&
                 fabs(fabs(dist[j]) - fabs(localdist[j])) < tol){
                dist[j] = std::max(dist[j], localdist[j]);
              }
              else{
                dist[j] = (fabs(dist[j]) < fabs(localdist[j])) ? dist[j] : localdist[j];
              }
            }
          }
        }
      }
      for (unsigned int j = 0; j < dist.size(); j++)
        box.setNodalValue(indices[j], dist[j]);
    
      if(box.getChildBox()) computeLevelset(gm, *box.getChildBox());
    }
    
    static void fillPointCloud(GEdge *ge, double sampling, std::vector<SPoint3> &points)
    {
      Range<double> t_bounds = ge->parBounds(0);
      double t_min = t_bounds.low();
      double t_max = t_bounds.high();
      double length = ge->length(t_min, t_max, 20);
      int N = length / sampling;
      for(int i = 0; i < N; i++) {
        double t = t_min + (double)i / (double)(N - 1) * (t_max - t_min);
        GPoint p = ge->point(t);
        points.push_back(SPoint3(p.x(), p.y(), p.z()));
      }
    }
    
    static int removeBadChildCells(cartesianBox<double> *parent)
    {
      cartesianBox<double> *child = parent->getChildBox();
      if(!child) return 0;
      int I = parent->getNxi();
      int J = parent->getNeta();
      int K = parent->getNzeta();
      for(int i = 0; i < I; i++)
        for(int j = 0; j < J; j++)
          for(int k = 0; k < K; k++){
            // remove children if they do not entirely fill parent, or if
            // there is no parent on the boundary
            int idx[8] = {child->getCellIndex(2 * i, 2 * j, 2 * k),
                          child->getCellIndex(2 * i, 2 * j, 2 * k + 1),
                          child->getCellIndex(2 * i, 2 * j + 1, 2 * k),
                          child->getCellIndex(2 * i, 2 * j + 1, 2 * k + 1),
                          child->getCellIndex(2 * i + 1, 2 * j, 2 * k),
                          child->getCellIndex(2 * i + 1, 2 * j, 2 * k + 1),
                          child->getCellIndex(2 * i + 1, 2 * j + 1, 2 * k),
                          child->getCellIndex(2 * i + 1, 2 * j + 1, 2 * k + 1)};
            bool exist[8], atLeastOne = false, butNotAll = false;
            for(int ii = 0; ii < 8; ii++){
              exist[ii] = child->activeCellExists(idx[ii]);
              if(exist[ii]) atLeastOne = true;
              if(!exist[ii]) butNotAll = true;
            }
            if(atLeastOne &&
               (butNotAll ||
                (i != 0 && !parent->activeCellExists(parent->getCellIndex(i - 1, j, k))) ||
                (i != I - 1 && !parent->activeCellExists(parent->getCellIndex(i + 1, j, k))) ||
                (j != 0 && !parent->activeCellExists(parent->getCellIndex(i, j - 1, k))) ||
                (j != J - 1 && !parent->activeCellExists(parent->getCellIndex(i, j + 1, k))) ||
                (k != 0 && !parent->activeCellExists(parent->getCellIndex(i, j, k - 1))) ||
                (k != K - 1 && !parent->activeCellExists(parent->getCellIndex(i, j, k + 1)))))
                for(int ii = 0; ii < 8; ii++) child->eraseActiveCell(idx[ii]);
          }
      return removeBadChildCells(child);
    }
    
    static void removeParentCellsWithChildren(cartesianBox<double> *box)
    {
      if(!box->getChildBox()) return;
      for(int i = 0; i < box->getNxi(); i++)
        for(int j = 0; j < box->getNeta(); j++)
          for(int k = 0; k < box->getNzeta(); k++){
            if(box->activeCellExists(box->getCellIndex(i, j, k))){
              cartesianBox<double> *parent = box, *child;
              int ii = i, jj = j, kk = k;
              while((child = parent->getChildBox())){
                ii *= 2; jj *= 2; kk *= 2;
                if(child->activeCellExists(child->getCellIndex(ii, jj, kk))){
                  box->eraseActiveCell(box->getCellIndex(i, j, k));
                  break;
                }
                parent = child;
              }
            }
          }
      removeParentCellsWithChildren(box->getChildBox());
    }
    
    static void removeOutsideCells(cartesianBox<double> *box)
    {
      if(!box) return;
      int nbErased = 0;
      for(cartesianBox<double>::cellIter it = box->activeCellsBegin();
          it != box->activeCellsEnd();){
        std::vector<double> vals;
        box->getNodalValuesForCell(*it, vals);
        double lsmax = *std::max_element(vals.begin(), vals.end());
        double lsmin = *std::min_element(vals.begin(), vals.end());
        double change_sign = lsmax * lsmin;
        double epsilon = 1.e-10;
        if(change_sign > 0 && lsmax < -epsilon) {
          box->eraseActiveCell(*(it++));
          nbErased++;
        }
        else ++it;
      }
      Msg::Info("  number of cells erased after filtering: %d", nbErased);
      removeOutsideCells(box->getChildBox());
    }
    
    int main(int argc,char *argv[])
    {
      if(argc < 6){
        printf("Usage: %s file lx ly lz rmax [levels=1] [refcs=1]\n", argv[0]);
        printf("where\n");
        printf("  'file' contains a CAD model\n");
        printf("  'lx', 'ly' and 'lz' are the sizes of the elements along the"
               " x-, y- and z-axis at the coarsest level\n");
        printf("  'rmax' is the radius of the largest sphere that can be inscribed"
               " in the structure\n");
        printf("  'levels' sets the number of levels in the grid\n");
        printf("  'refcs' selects if curved surfaces should be refined\n");
        return -1;
      }
    
      GmshInitialize();
      GmshSetOption("General", "Terminal", 1.);
      GmshMergeFile(argv[1]);
      double lx = atof(argv[2]), ly = atof(argv[3]), lz = atof(argv[4]);
      double rmax = atof(argv[5]);
      int levels = (argc > 6) ? atof(argv[6]) : 1;
      int refineCurvedSurfaces = (argc > 7) ? atof(argv[7]) : 1;
    
      // minimum distance between points in the cloud at the coarsest
      // level
      double sampling = std::min(rmax, std::min(lx, std::min(ly, lz)));
    
      // radius of the "tube" created around parts to refine at the
      // coarsest level
      double rtube = std::max(lx, std::max(ly, lz)) * 2.;
    
      GModel *gm = GModel::current();
    
      std::vector<SPoint3> points;
      Msg::Info("Filling coarse point cloud on surfaces");
      for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++)
        (*fit)->fillPointCloud(sampling, &points);
      Msg::Info("  %d points in the surface cloud", (int)points.size());
    
      std::vector<SPoint3> refinePoints;
      if(levels > 1){
        double s = sampling / pow(2., levels - 1);
        Msg::Info("Filling refined point cloud on curves and curved surfaces");
        for (GModel::eiter eit = gm->firstEdge(); eit != gm->lastEdge(); eit++)
          fillPointCloud(*eit, s, refinePoints);
    
        // FIXME: refine this by computing e.g. "mean" curvature
        if(refineCurvedSurfaces){
          for (GModel::fiter fit = gm->firstFace(); fit != gm->lastFace(); fit++)
            if((*fit)->geomType() != GEntity::Plane)
              (*fit)->fillPointCloud(2 * s, &refinePoints);
        }
        Msg::Info("  %d points in the refined cloud", (int)refinePoints.size());
      }
    
      SBoundingBox3d bb;
      for(unsigned int i = 0; i < points.size(); i++) bb += points[i];
      for(unsigned int i = 0; i < refinePoints.size(); i++) bb += refinePoints[i];
      bb.scale(1.21, 1.21, 1.21);
      SVector3 range = bb.max() - bb.min();
      int NX = range.x() / lx;
      int NY = range.y() / ly;
      int NZ = range.z() / lz;
      if(NX < 2) NX = 2;
      if(NY < 2) NY = 2;
      if(NZ < 2) NZ = 2;
    
      Msg::Info("  bounding box min: %g %g %g -- max: %g %g %g",
                bb.min().x(), bb.min().y(), bb.min().z(),
                bb.max().x(), bb.max().y(), bb.max().z());
      Msg::Info("  Nx=%d Ny=%d Nz=%d", NX, NY, NZ);
    
      cartesianBox<double> box(bb.min().x(), bb.min().y(), bb.min().z(),
                               SVector3(range.x(), 0, 0),
                               SVector3(0, range.y(), 0),
                               SVector3(0, 0, range.z()),
                               NX, NY, NZ, levels);
    
      Msg::Info("Inserting active cells in the cartesian grid");
      Msg::Info("  level %d", box.getLevel());
      for (unsigned int i = 0; i < points.size(); i++)
        insertActiveCells(points[i].x(), points[i].y(), points[i].z(), rmax, box);
    
      cartesianBox<double> *parent = &box, *child;
      while((child = parent->getChildBox())){
        Msg::Info("  level %d", child->getLevel());
        for(unsigned int i = 0; i < refinePoints.size(); i++)
          insertActiveCells(refinePoints[i].x(), refinePoints[i].y(), refinePoints[i].z(),
                            rtube / pow(2., (levels - child->getLevel())), *child);
        parent = child;
      }
    
      // remove child cells that do not entirely fill parent cell or for
      // which there is no parent neighbor; then remove parent cells that
      // have children
      Msg::Info("Removing cells to match X-FEM mesh topology constraints");
      removeBadChildCells(&box);
      removeParentCellsWithChildren(&box);
    
      // we generate duplicate nodes at this point so we can easily access
      // cell values at each level; we will clean up by renumbering after
      // filtering
      Msg::Info("Initializing nodal values in the cartesian grid");
      box.createNodalValues();
    
      Msg::Info("Computing levelset on the cartesian grid");
      computeLevelset(gm, box);
    
      Msg::Info("Removing cells outside the structure");
      removeOutsideCells(&box);
    
      Msg::Info("Renumbering mesh vertices across levels");
      box.renumberNodes();
    
      bool decomposeInSimplex = false;
      box.writeMSH("yeah.msh", decomposeInSimplex);
    
      Msg::Info("Done!");
      GmshFinalize();
    }