ibgext.c

有限元学习研究用源代码(老外的),供科研人员参考

#include "ibg.h"
#include "ibgg.h"
#include "ibglib.h"
#include "ibgerror.h"
#include "ibgext.h"
#include <stdio.h>
#include <math.h>

#define MaxFunctions 50
#define MaxValues 50

struct ibgExtData0{
    ibgFloat  X[ibgDIM][MaxValues];
    ibgInt    Num[ibgDIM];
    ibgInt    Reg[ibgDIM];
    ibgFloat  Min;
    ibgFloat  Dmin[ibgDIM];
    ibgFloat  FaceNormal;
    ibgFloat  LineNormal;
    ibgFloat  FaceTangential;
    ibgFloat  LineTangential;
    ibgInt    Type[MaxFunctions];
    ibgInt    Dir[MaxFunctions];
    ibgInt    Region[MaxFunctions];
    ibgFloat  Value[MaxFunctions];
    ibgFloat  Pos[MaxFunctions][6];
    ibgFloat  Fac[MaxFunctions][6];
    int       Number;
};

ibgInt  ibgExtReg(ibgExtData *r, ibgInt d)
{return r->Reg[d];}
ibgInt  ibgExtNum(ibgExtData *r, ibgInt d)
{return r->Num[d];}
ibgFloat* ibgExtX(ibgExtData *r, ibgInt d)
{return &(r->X[d][0]);}

static ibgExtData *ibgExt;
	      
char *names[]={"minX","minY","minZ","min",
	       "regularX","regularY","regularZ","regular",
	       "tangentialToFace","normalToFace",
	       "tangentialToLine","normalToLine",
	       "refineX","refineY","refineZ","refineU","refineV","refineW",
	       "refine",
	       "x","y","z",
	       ""};
int   scalar[]={1,1,1,1,
		1,1,1,1,
		1,1,
		1,1,
		0,
		0,0,0,0,0,0,
		0,0,0,
		2};
/* for every scalar and nonscalar name there is one case below.
   be careful with the related numbers if you change this.
   For example, include a name beginning with x before "x", but don't
   forget to move the value for "x" in the related case statement.
*/
#define buflen 200
static char buffer[buflen];
static const ibgFloat verybig = 1.e36;
static const ibgFloat verysmall = 1.e-36;
static const ibgFloat nonsense = 1.1e36;
static int b; 
static FILE *f;
	      
static int nextchar()
{char* rc;      
  beg:	      
    while((buffer[b] == ' ')||(buffer[b]=='\t')) b++;
    if( (buffer[b] == '!' ) || (buffer[b] == '?' )   ||
        (buffer[b] == '%' ) || (buffer[b] == '#' )   ||
        (buffer[b] == '\n') || (buffer[b] == '\0')   ||
       ((buffer[b] == '-' ) && (buffer[b+1] == '-')) ||
       ((buffer[b] == '/' ) && (buffer[b+1] == '/')) ){
        rc=fgets(buffer,buflen,f); b=0;
        if(rc==NULL){b=0; buffer[b] = '\0'; return 0;}
	goto beg;
    }	      
    return buffer[b];
}	      
	      
static ibgFloat nextReal()
{ibgFloat x; int rc;
 rc = sscanf(&buffer[b],"%f",&x);
 if(rc!=1) return nonsense;
  ignore:     
 if((buffer[b]=='+')||(buffer[b]=='-')){b++;}
 while((buffer[b]>='0')&&(buffer[b]<='9')) b++;
 if((buffer[b]=='.')){b++;}
 while((buffer[b]>='0')&&(buffer[b]<='9')) b++;
 if((buffer[b]=='e')){b++;}
 if((buffer[b]=='+')||(buffer[b]=='-')){b++;}
 while((buffer[b]>='0')&&(buffer[b]<='9')) b++;
 nextchar();  
 return x;    
}	      
	      
ibgExtData *ibgExtInit(char* filename)
{int i,j,k,n,t,d,xi;
 char* rc;
 ibgFloat x, xx[MaxValues];   
 ibgExt = (ibgExtData *)malloc(sizeof(ibgExtData));
 ibgExt->Number = 0;
 ibgExt->Min     = verysmall;
 ibgExt->Dmin[0] = verysmall;
 ibgExt->Dmin[1] = verysmall;
 ibgExt->Dmin[2] = verysmall;
 ibgExt->Num[0] = 0; 
 ibgExt->Num[1] = 0; 
 ibgExt->Num[2] = 0; 
 ibgExt->Reg[0] = 0; 
 ibgExt->Reg[1] = 0; 
 ibgExt->Reg[2] = 0; 
 ibgExt->X[0][0] = 0;
 ibgExt->X[1][0] = 0;
 ibgExt->X[2][0] = 0;
 ibgExt->FaceNormal = verybig;
 ibgExt->LineNormal = verybig;
 ibgExt->FaceTangential = verybig;
 ibgExt->LineTangential = verybig;
 f = fopen(filename,"r"); 
 if(f==NULL) {
     ibgmessage(ibgMEOpenFile);
     ibgprintf("%s",filename);
     ibgmsgreset;
 }
 b=0; buffer[0]='\0';
 while(nextchar()){
     for(i=0;names[i][0]!='\0';i++){
         for(j=0;;j++){
	     if(names[i][j]=='\0') goto found;
             if(names[i][j]!=buffer[b+j]) break;
         }    
     }	      
     goto incformat;
   found:     
     b += j;  
     if(scalar[i]){ 
         nextchar();
	 if(buffer[b]!='=') goto incformat;
	 b++; nextchar();
	 x = nextReal(); xi = (int)(x+0.1);
	 if(x==nonsense) goto incformat;
	 switch(i){
	 case 0:          /* minX */
             ibgExt->Dmin[0]=x;
	     break;
	 case 1:          /* minY */
	     ibgExt->Dmin[1]=x;
	     break;
	 case 2:          /* minZ */
	     ibgExt->Dmin[2]=x;
	     break;
	 case 3:          /* min */
	     ibgExt->Dmin[0]=ibgExt->Dmin[1]=ibgExt->Dmin[2]=x;
	     break;
	 case 4:          /* regularX */
             ibgExt->Reg[0]=xi;
	     break;
	 case 5:          /* regularY */
	     ibgExt->Reg[1]=xi;
	     break;
	 case 6:          /* regularZ */
	     ibgExt->Reg[2]=xi;
	     break;
	 case 7:          /* regular */
	     ibgExt->Reg[0]=ibgExt->Reg[1]=ibgExt->Reg[2]=xi;
	     break;
	 case 8:          /* tangentialToFace */
	     ibgExt->FaceTangential=x;
	     break;
	 case 9:          /* normalToFace */
	     ibgExt->FaceNormal=x;
	     break;
	 case 10:          /* tangentialToLine */
	     ibgExt->LineTangential=x;
	     break;
	 case 11:          /* normalToLine */
	     ibgExt->LineNormal=x;
	     break;
	 }
     }else{   
	 nextchar();
	 for(k=0;k<MaxValues;k++) xx[k] = 0;
         if(buffer[b]!='(') goto incformat;
	 b++; nextchar();
	 k=0; 
	 while(1){
	     while(buffer[b]==','){
                 k++; b++;
                 nextchar();
	     }
	     if(buffer[b]==')') {b++;break;}
	     xx[k] = nextReal();
	     if(xx[k]==nonsense) {goto incformat;}
             nextchar();
	 }    
	 if((n=ibgExt->Number)>=MaxFunctions) goto incformat;
	 switch(i){
	 case 12:          /* refineX */
	     t = 1;
	     d = 0;
	     break;
	 case 13:          /* refineY */
	     t = 1;
	     d = 1;
	     break;
	 case 14:          /* refineZ */
	     t = 1;
	     d = 2;
	     break;
	 case 15:          /* refineU */
	     t = 1;
	     d = 3;
	     break;
	 case 16:          /* refineV */
	     t = 1;
	     d = 4;
	     break;
	 case 17:          /* refineW */
	     t = 1;
	     d = 5;
	     break;
	 case 18:          /* refine */
	     t = 0;
	     break;
	 case 19:          /* x */
	     t = 2;
	     d = 0;
	     break;
	 case 20:          /* y */
	     t = 2;
	     d = 1;
	     break;
	 case 21:          /* z */
	     t = 2;
	     d = 2;
	     break;
	 }
	 switch(t){
	 case 0:
	 case 1:
	     ibgExt->Type[n]   = t;
	     ibgExt->Dir[n]    = d;
	     ibgExt->Value[n]  = xx[0]*xx[0];
	     ibgExt->Region[n] = (int) (xx[1]+0.1);
	     ibgExt->Pos[n][0] = xx[2];
	     ibgExt->Fac[n][0] = xx[3];
	     ibgExt->Pos[n][1] = xx[4];
	     ibgExt->Fac[n][1] = xx[5];
	     ibgExt->Pos[n][2] = xx[6];
	     ibgExt->Fac[n][2] = xx[7];
	     ibgExt->Pos[n][3] = xx[8];
	     ibgExt->Fac[n][3] = xx[9];
	     ibgExt->Pos[n][4] = xx[10];
	     ibgExt->Fac[n][4] = xx[11];
	     ibgExt->Pos[n][5] = xx[12];
	     ibgExt->Fac[n][5] = xx[13];
	     if(ibgExt->Value[n]<ibgExt->Min*ibgExt->Min) 
		 ibgExt->Value[n]=ibgExt->Min*ibgExt->Min;
	     if(ibgExt->Region[n]<0) ibgExt->Region[n]=0;
	     if(ibgExt->Region[n]>=ibgSMAX) ibgExt->Region[n]=0;
	     ibgExt->Number++;
	     break;
	 case 2:
	     ibgExt->Num[d] = k+1; x = -verybig;
	     for(j=0;j<=k;j++){
		 if(x>=xx[j]){ibgExt->Num[d] = j; break;}
		 ibgExt->X[d][j] = x = xx[j];
	     }
	     if( ibgExt->Num[d] < 0){
		 ibgExt->Num[d] = 0; ibgExt->X[d][j] = 0;
	     }
	 }
     }
     nextchar();
     if((buffer[b]!=',')&&(buffer[b]!=';')) goto incformat;
     b++;
 }	      
 goto afterread;
  incformat:  
 ibgmessage(ibgMEWrongFormat);
 ibgprintf("%s\n",filename);
 ibgprintf("%s",buffer);
 for(i=0;i<b;i++) ibgprintf("-");
 ibgprintf(">");
 ibgmsgreset;
  afterread:
 fclose(f);
 ibggDmin[0] = ibgExt->Dmin[0];
 ibggDmin[1] = ibgExt->Dmin[1];
 ibggDmin[2] = ibgExt->Dmin[2];
 ibggDefaultFaceNormal = ibgExt->FaceNormal;
 ibggDefaultLineNormal = ibgExt->LineNormal;
 ibggDefaultFaceTangential = ibgExt->FaceTangential;
 ibggDefaultLineTangential = ibgExt->LineTangential;
 return ibgExt;
}


void ibgExtRegion(ibgPtObject This, ibgPoint *point,
        		   	      	ibgFloat *length)
{ibgExtData *r = (ibgExtData *) This;
 int i; ibgFloat rr,x,y,z;
 *length = verybig;
 for(i=0;i<r->Number;i++){
     if(r->Type[i]!=0) continue;
     if((r->Region[i]) &&(r->Region[i] != ibgpSegment(*point)))continue;
     x = ibgpX(*point)[0] - r->Pos[i][0];
     y = ibgpX(*point)[1] - r->Pos[i][1];
     z = ibgpX(*point)[2] - r->Pos[i][2];
     x *= r->Fac[i][0];
     y *= r->Fac[i][1];
     z *= r->Fac[i][2];
     rr = (ibgFloat)sqrt((double)(r->Value[i] + x*x + y*y + z*z));
     if(rr<*length) *length=rr;
 }
}

int ibgExtEdge(ibgPtObject This, ibgPoint *x1, ibgPoint *x2)
{ibgExtData *r = (ibgExtData *) This;
 int i; ibgFloat rr,dd,x,y,z;
 ibgFloat x0=(ibgpX(*x1)[0]+ibgpX(*x2)[0])/2;
 ibgFloat y0=(ibgpX(*x1)[1]+ibgpX(*x2)[1])/2;
 ibgFloat z0=(ibgpX(*x1)[2]+ibgpX(*x2)[2])/2;
 for(i=0;i<r->Number;i++){
     if     (r->Type[i]!=1) continue;
     if((r->Region[i]) &&
	((r->Region[i] != ibgpSegment(*x1)) &&
	 (r->Region[i] != ibgpSegment(*x2)))  ) continue;
     dd = ibgpX(*x1)[r->Dir[i]]-ibgpX(*x2)[r->Dir[i]];
     if(dd<0) dd = -dd;
     x = x0 - r->Pos[i][0];
     y = y0 - r->Pos[i][1];
     z = z0 - r->Pos[i][2];
     x *= r->Fac[i][0];
     y *= r->Fac[i][1];
     z *= r->Fac[i][2];
     rr = (ibgFloat)sqrt((double)(r->Value[i] + x*x + y*y + z*z));
     if(rr<dd) return ibgTrue;
 }
 return ibgFalse;
}