Skip to content
Snippets Groups Projects
Select Git revision
  • gmsh_4_12_2
  • 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_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
40 results

x4.jl

Blame
  • OctreeInternals.cpp 13.38 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 <list>
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include "GmshMessage.h"
    #include "OctreeInternals.h"
    
    int initializeOctantBuckets(double *_orig, double *_size, int _maxElem,
                                octantBucket **buckets_head,
                                globalInfo **globalPara)
    // Initialize the buckets
    // Given by user: orig and size -- information about the domain
    //                maxElem -- maximum number of elements per bucket
    // Return: buckets -- pointer to the begin of buckets
    //         globalPara -- some info about the buckets
    // At last, 1 will be returned if succeed, otherwise, return 0
    {
      int i, j, k, tmp1;
      int p = 1;
      int initial_buckets_num; // which is a number of 8^p form for integer p
      double tmp[3], error[3];
      octantBucket *buckets;
    
      for(i = 0; i < 3; i++) error[i] = _size[i] * 0.01;
    
      initial_buckets_num = (int)pow(8., p); // it is actually 8
    
      (*globalPara) = new globalInfo;
      (*globalPara)->maxPrecision = 1;
      (*globalPara)->maxElements = _maxElem;
      (*globalPara)->ptrToPrevElement = NULL;
    
      for(i = 0; i < 3; i++) {
        (*globalPara)->origin[i] = _orig[i];
        (*globalPara)->size[i] = _size[i];
      }
    
      (*globalPara)->numBuckets = initial_buckets_num;
      *buckets_head = new octantBucket;
      if(!(*buckets_head)) {
        Msg::Error("initializeOctantBuckets could not allocate enough space");
        return (0);
      } // if could not allocate buckets
    
      buckets = new octantBucket[8];
      if(!buckets) {
        Msg::Error("initializeOctantBuckets could not allocate enough space");
        return (0);
      }
    
      (*buckets_head)->next = buckets;
      (*buckets_head)->parent = NULL;
      (*buckets_head)->numElements = 0;
      (*buckets_head)->lhead = NULL;
      (*buckets_head)->precision = 0;
      for(i = 0; i < 3; i++) {
        (*buckets_head)->minPt[i] = _orig[i] - error[i];
        (*buckets_head)->maxPt[i] = _size[i] + _orig[i] + error[i];
      }
    
      for(i = 0; i < (*globalPara)->numBuckets; i++) {
        buckets[i].numElements = 0;
        buckets[i].lhead = NULL;
        buckets[i].next = NULL;
        buckets[i].parent = *buckets_head;
        buckets[i].precision = 1;
      }
    
      tmp1 = (int)(pow(2., p));
      for(i = 0; i < 3; i++) {
        tmp[i] = (double)(_size[i] + 2 * error[i]) / tmp1;
      }
    
      for(k = 0; k < tmp1; k++) {
        for(j = 0; j < tmp1; j++) {
          for(i = 0; i < tmp1; i++) {
            buckets[i + j * tmp1 + k * tmp1 * tmp1].minPt[0] =
              (*buckets_head)->minPt[0] + tmp[0] * i;
            buckets[i + j * tmp1 + k * tmp1 * tmp1].minPt[1] =
              (*buckets_head)->minPt[1] + tmp[1] * j;
            buckets[i + j * tmp1 + k * tmp1 * tmp1].minPt[2] =
              (*buckets_head)->minPt[2] + tmp[2] * k;
            buckets[i + j * tmp1 + k * tmp1 * tmp1].maxPt[0] =
              (*buckets_head)->minPt[0] + tmp[0] * (i + 1);
            buckets[i + j * tmp1 + k * tmp1 * tmp1].maxPt[1] =
              (*buckets_head)->minPt[1] + tmp[1] * (j + 1);
            buckets[i + j * tmp1 + k * tmp1 * tmp1].maxPt[2] =
              (*buckets_head)->minPt[2] + tmp[2] * (k + 1);
          }
        }
      }
    
    #if 0
      for (i = 0; i < 8; i++) {
        printf(" bucket %d : min[0]=%f, min[1]=%f, min[2]=%f, max[0]= %f, max[1]=%f, max[3]=%f\n",
               i,buckets[i].minPt[0],buckets[i].minPt[1],
               buckets[i].minPt[2], buckets[i].maxPt[0], buckets[i].maxPt[1],
               buckets[i].maxPt[2]);
        printf("bucket elements link list: bucket->lhead = %d\n", buckets[i].lhead);
      }
    #endif
    
      return (1);
    }
    
    int addElement2Bucket(octantBucket *_bucket, void *_element, double *_minBB,
                          double *_maxBB, double *_ele_centroid,
                          globalInfo *_globalPara)
    // Add another element to the octant bucket's list.
    // If the bucket contains too many elements after adding this element,
    // refine this bucket and reallocate the elements of this bucket
    // Given:- the octant bucket, - the element
    //       - the element's minimum and maximum x,y,z
    //       - the element's centroid,  - global information
    // Check if element has already been added - if not, return 1
    // for successfully adding, otherwise return -1
    {
      int i, flag = 1;
      ELink ptr1, ptr2;
      octantBucket *ptrBucket;
    
      // check for duplicates
      if(checkElementInBucket(_bucket, _element) == 1) return -1;
    
      // printf("\n addToBucket...\n");
      // ptr1 = (ELink) malloc(sizeof(Elem));
      ptr1 = new Elem;
      (_globalPara->listAllElements).push_back(_element);
    
      ptr1->next = _bucket->lhead;
      ptr1->region = _element;
    
      for(i = 0; i < 3; i++) {
        ptr1->minPt[i] = _minBB[i];
        ptr1->maxPt[i] = _maxBB[i];
        ptr1->centroid[i] = _ele_centroid[i];
        // printf(" %7.2f->%-7.2f",ptr1->minPt[i], ptr1->maxPt[i]);
      }
    
      _bucket->lhead = ptr1;
      (_bucket->numElements)++;
    
    #if 0
      printf("bucket element list: bucket->lhead = %d", _bucket->lhead);
      printf(" numElements = %d\n",_bucket->numElements);
      printf("the element is add to this bucket: (%f, %f, %f) to (%f, %f, %f)\n",
             _bucket->minPt[0],_bucket->minPt[1], _bucket->minPt[2], _bucket->maxPt[0],
             _bucket->maxPt[1], _bucket->maxPt[2]);
    #endif
    
      // check whether the number of elements in the bucket > maxElements
      // if true, refine the bucket and reallocate the elements
      while(flag == 1) {
        flag = 0;
        if(_bucket->numElements > _globalPara->maxElements) {
          // printf(" going to subdivide\n");
    
          subdivideOctantBucket(_bucket, _globalPara);
    
          // printf("finish subdivede \n");
    
          ptr1 = _bucket->lhead;
          while(ptr1 != NULL) {
            ptrBucket = findElementBucket(_bucket, ptr1->centroid);
            ptr2 = ptr1;
            ptr1 = ptr1->next;
            if(ptrBucket == NULL) {
              Msg::Error("Wrong , ptrBucket = NULL. A bug here!");
              return 0;
            }
            ptr2->next = ptrBucket->lhead;
            ptrBucket->lhead = ptr2;
            (ptrBucket->numElements)++;
            if(ptrBucket->numElements > _globalPara->maxElements) {
              flag = 1;
              _bucket->lhead = NULL;
              _bucket = ptrBucket;
            }
          }
          if(flag == 0) _bucket->lhead = NULL;
        }
      }
      return 1;
    }
    
    int checkElementInBucket(octantBucket *_bucket, void *_element)
    // Given an elememt and an octant bucket, check if the element
    // exists in the bucket's element list. return 1 if already exits,
    // otherwise, return 0
    {
      ELink ptr;
      for(ptr = _bucket->lhead; ptr != NULL; ptr = ptr->next) {
        // changed ****, compare the objected pointed by the void *.
        if(ptr->region == _element) return 1;
      }
      return 0;
    }
    
    octantBucket *findElementBucket(octantBucket *_buckets_head, double *_pt)
    // Find the leaf bucket which contains the point _pt
    // given parameter: _buckets --- the point to buckets head
    //                  _pt --- the point to find
    // Return the pointer to the bucket contains the point
    // if fail, return NULL
    {
      int i, j;
      int num = 8;
      octantBucket *prevbucket = NULL;
      octantBucket *tmpbucket = _buckets_head->next;
    
      while(tmpbucket != NULL) {
        for(i = 0; i < num; i++) {
          for(j = 0; j < 3; j++) {
            if(tmpbucket[i].minPt[j] > _pt[j] || tmpbucket[i].maxPt[j] < _pt[j])
              break;
          }
          if(j == 3) {
            prevbucket = tmpbucket + i;
            tmpbucket = tmpbucket[i].next;
            break;
          }
        } // for loop i
        if(i == num) {
          // Msg::Error("No bucket contains the given point!");
          return NULL;
        }
      } // for while loop
      return prevbucket;
    }
    
    int subdivideOctantBucket(octantBucket *_bucket, globalInfo *_globalPara)
    // To many elements are in this octant bucket, so try to refine
    // Returns 1 for success, 0 for failure (no memory left).
    {
      int i, j, k, tmp1;
      int numBuck = 8;
      double tmp[3];
    
      _bucket->next = new octantBucket[8];
      // _bucket->next  = (octantBucket *) calloc(numBuck,sizeof(octantBucket));
    
      if(!_bucket->next) {
        Msg::Error("subdivideOctantBucket could not allocate enough space");
        return 0;
      }
    
      _globalPara->numBuckets += 8;
      if(_bucket->precision == _globalPara->maxPrecision)
        _globalPara->maxPrecision++;
      for(i = 0; i < numBuck; i++) {
        (_bucket->next[i]).numElements = 0;
        (_bucket->next[i]).lhead = NULL;
        (_bucket->next[i]).next = NULL;
        (_bucket->next[i]).parent = _bucket;
        (_bucket->next[i]).precision = _bucket->precision + 1;
      }
    
      tmp1 = 2;
      for(i = 0; i < 3; i++) {
        tmp[i] = ((double)(_bucket->maxPt[i] - _bucket->minPt[i])) / tmp1;
      }
    
      for(k = 0; k < tmp1; k++) {
        for(j = 0; j < tmp1; j++) {
          for(i = 0; i < tmp1; i++) {
            _bucket->next[i + j * tmp1 + k * tmp1 * tmp1].minPt[0] =
              _bucket->minPt[0] + tmp[0] * i;
            _bucket->next[i + j * tmp1 + k * tmp1 * tmp1].minPt[1] =
              _bucket->minPt[1] + tmp[1] * j;
            _bucket->next[i + j * tmp1 + k * tmp1 * tmp1].minPt[2] =
              _bucket->minPt[2] + tmp[2] * k;
            _bucket->next[i + j * tmp1 + k * tmp1 * tmp1].maxPt[0] =
              _bucket->minPt[0] + tmp[0] * (i + 1);
            _bucket->next[i + j * tmp1 + k * tmp1 * tmp1].maxPt[1] =
              _bucket->minPt[1] + tmp[1] * (j + 1);
            _bucket->next[i + j * tmp1 + k * tmp1 * tmp1].maxPt[2] =
              _bucket->minPt[2] + tmp[2] * (k + 1);
          }
        }
      }
    
      return 1;
    }
    
    void *searchElement(octantBucket *_buckets_head, double *_pt,
                        globalInfo *_globalPara, BBFunction BBElement,
                        InEleFunction xyzInElement)
    {
      int flag;
      octantBucket *ptrBucket;
      ELink ptr1;
      std::vector<void *>::iterator iter;
      void *ptrToEle = _globalPara->ptrToPrevElement;
    
      if(ptrToEle) {
        flag = xyzInElementBB(_pt, ptrToEle, BBElement);
        if(flag == 1) flag = xyzInElement(ptrToEle, _pt);
        if(flag == 1) return ptrToEle;
      }
    
      ptrBucket = findElementBucket(_buckets_head, _pt);
      if(ptrBucket == NULL) {
        // this is not an error
        Msg::Debug("The point is not in the domain");
        return NULL;
      }
    
      ptr1 = ptrBucket->lhead;
    
    #if 0
      printf("point %lf %lf %lf has been found in bucket %lf %lf %fl -> %lf %lf %lf  %p\n",
             _pt[0],_pt[1],_pt[2], ptrBucket->minPt[0],ptrBucket->minPt[1],ptrBucket->minPt[2],
             ptrBucket->maxPt[0],ptrBucket->maxPt[1],ptrBucket->maxPt[2], ptr1);
      if (ptr1 == NULL) {
        printf("empty element list for centroid list!?\n, possible!");
      }
    #endif
    
      while(ptr1 != NULL) {
        flag = xyzInElementBB(_pt, ptr1->region, BBElement);
        if(flag == 1) flag = xyzInElement(ptr1->region, _pt);
        if(flag == 1) {
          _globalPara->ptrToPrevElement = ptr1->region;
          return ptr1->region;
        }
        ptr1 = ptr1->next;
      }
    
      for(iter = (ptrBucket->listBB).begin(); iter != (ptrBucket->listBB).end();
          iter++) {
        flag = xyzInElementBB(_pt, *iter, BBElement);
        if(flag == 1) flag = xyzInElement(*iter, _pt);
        if(flag == 1) {
          _globalPara->ptrToPrevElement = *iter;
          return *iter;
        }
      }
    
      // printf("This point is not found in all elements! It is not in the domain
      // \n");
      return NULL;
    }
    
    int xyzInElementBB(double *_xyz, void *_region, BBFunction _bbElement)
    // Check if xyz is in the region's bounding box, return 1 if true, 0 otherwise
    // BBElement is the function given by user to find the bounding box
    {
      int i;
      double minPt[3]; // corner with smallest x,y,z coords
      double maxPt[3]; // corner with largest x,y,z coords
    
      (*_bbElement)(_region, minPt, maxPt);
    
      for(i = 0; i < 3; i++) {
        if(_xyz[i] > maxPt[i] || _xyz[i] < minPt[i]) return 0;
      } // for ith direction
    
      return 1;
    }
    
    void insertOneBB(void *_region, double *_minPt, double *_maxPt,
                     octantBucket *_bucket)
    {
      int i;
      ELink ptr;
      for(i = 0; i < 3; i++) {
        if(_bucket->minPt[i] > _maxPt[i] || _bucket->maxPt[i] < _minPt[i]) return;
      }
      if(_bucket->next == NULL) {
        ptr = _bucket->lhead;
        while(ptr != NULL) {
          if(ptr->region == _region) return;
          ptr = ptr->next;
        }
    
        //_bucket->listBB.insert(_bucket->listBB.end(),_region);
        _bucket->listBB.push_back(_region);
        return;
      }
    
      for(i = 0; i < 8; i++)
        insertOneBB(_region, _minPt, _maxPt, _bucket->next + i);
      return;
    }
    
    void *searchAllElements(octantBucket *_buckets_head, double *_pt,
                            globalInfo *_globalPara, BBFunction BBElement,
                            InEleFunction xyzInElement,
                            std::vector<void *> *_elements)
    {
      int flag, flag1;
      octantBucket *ptrBucket;
      std::vector<void *>::iterator iter;
    
      ptrBucket = findElementBucket(_buckets_head, _pt);
      if(ptrBucket == NULL) {
        Msg::Debug("The point is not in the domain");
        return NULL;
      }
    
    #if 0
      printf("point %lf %lf %lf has been found in bucket %lf %lf %fl -> %lf %lf %lf  %p\n",
             _pt[0],_pt[1],_pt[2], ptrBucket->minPt[0],ptrBucket->minPt[1],ptrBucket->minPt[2],
             ptrBucket->maxPt[0],ptrBucket->maxPt[1],ptrBucket->maxPt[2], ptr1);
    
      if (ptr1 == NULL) {
        printf("empty element list for centroid list!?\n, possible!");
      }
    #endif
    
      flag1 = 0;
      ELink ptr1 = ptrBucket->lhead;
      while(ptr1 != NULL) {
        flag = xyzInElementBB(_pt, ptr1->region, BBElement);
        if(flag == 1) flag = xyzInElement(ptr1->region, _pt);
        if(flag == 1) {
          _elements->push_back(ptr1->region);
          flag1 = 1;
        }
        ptr1 = ptr1->next;
      }
    
      for(iter = (ptrBucket->listBB).begin(); iter != (ptrBucket->listBB).end();
          iter++) {
        flag = xyzInElementBB(_pt, *iter, BBElement);
        if(flag == 1) flag = xyzInElement(*iter, _pt);
        if(flag == 1) {
          _elements->push_back(*iter);
          flag1 = 1;
        }
      }
    
      if(flag1) return (void *)(_elements);
    
      // Msg::Warning("This point is not found in any element! It is not in the
      // domain");
      return NULL;
    }