From 1d3b33037381455c71bc3ad21f5273256c4683a6 Mon Sep 17 00:00:00 2001 From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be> Date: Fri, 23 Apr 2004 17:02:13 +0000 Subject: [PATCH] *** empty log message *** --- Plugin/o_internals.cpp | 441 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 441 insertions(+) create mode 100644 Plugin/o_internals.cpp diff --git a/Plugin/o_internals.cpp b/Plugin/o_internals.cpp new file mode 100644 index 0000000000..0e94991bab --- /dev/null +++ b/Plugin/o_internals.cpp @@ -0,0 +1,441 @@ +#include <iostream> +#include <stdio.h> +#include <stdlib.h> +#include <math.h> +#include <list> + +using std::list; +#include "o_internals.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; //(globalInfo *)malloc(sizeof(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 = (octantBucket *)malloc(sizeof(octantBucket)); + *buckets_head = new octantBucket; + if (!(*buckets_head)) { + printf("Error, initializeOctantBuckets could not allocate enough space\n"); + return (0); + } /* if could not allocate buckets */ + +// buckets = (octantBucket *)calloc((*globalPara)->numBuckets,sizeof(octantBucket)); + buckets = new octantBucket[8]; + if (!buckets) { + printf("Error, initializeOctantBuckets could not allocate enough space\n"); + return (0); + } /* if could not allocate buckets */ + + (*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); + } + } + } + /* for test + 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); + } + end for test +*/ + + 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; /* flag to identify whether refine is finished or not */ + 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)++; + + /* + 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] ); + */ + + /* 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) + printf("Wrong , ptrBucket = NULL. A bug here!\n"); + 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) { + // printf("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) { + fprintf(stderr,"Error, subdivideOctantBucket could not allocate enough space\n"); + return 0; + } /* if could not allocate buckets */ + + _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; + list<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) { + // printf("Error! the point is not in the domain.\n"); + return NULL; + } + + ptr1 = ptrBucket->lhead; + /* + + 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!"); + } + */ + + 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 direciton */ + 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); + 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, list<void *> *_elements) +{ + int flag, flag1; + octantBucket *ptrBucket; + list<void *>::iterator iter; + + ptrBucket = findElementBucket(_buckets_head, _pt); + if(ptrBucket == NULL) { + printf("Error! the point is not in the domain.\n"); + return NULL; + } + + /* + 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!"); + } + */ + + flag1 = 0; + for(iter = (ptrBucket->listBB).begin(); iter != + (ptrBucket->listBB).end(); iter++) + { + printf("Enter 1 \n"); + 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); + + printf("This point is not found in all elements! It is not in the domain \n"); + return NULL; +} -- GitLab