00001 #include "meshtransfer.h" 00002 00003 #include "global.h" 00004 #include "globmat.h" 00005 #include "intpoints.h" 00006 #include "loadcase.h" 00007 #include "element.h" 00008 #include "plelemlt.h" 00009 #include "plelemqt.h" 00010 #include "plelemlq.h" 00011 #include "plelemqq.h" 00012 #include "nssolver.h" 00013 #include "stochdriver.h" 00014 #include "mefelinit.h" 00015 00016 #include "gtopology.h" 00017 #include "gelement.h" 00018 #include "gnode.h" 00019 #include "gmatrix.h" 00020 #include "mathem.h" 00021 #include "problem.h" 00022 //#include "least_square.h" 00023 00024 #include <time.h> 00025 #include <math.h> 00026 #include <stdlib.h> 00027 00028 //////////////////// /* termitovo */ //////////////////////////////////// 00029 00030 // /* 00031 // vse je delano pro 2D planestress prvky s ndofn=2 00032 // pro jine to neni kontrolovano !!!!!!!!!!!!!!!!!! 00033 // prvek plelemqq je predpokladan s rovnymi hranami(subparametricky), neb v fci 'nc_quad_4_2d' je zatim odkaz jen na 'nc_lin_4_2d' 00034 // */ 00035 // /* 00036 // ze site na sit se prenaseji posuny a prvnich ncompstr hodnot z pole other , ne prvnich ncompother neb gamma a hardening nepotrebuju 00037 // */ 00038 // /* 00039 // fce nezavisle na jine globalni promenne: 00040 // give_ndofn ok 00041 // give_dof ok 00042 // give_nip ok 00043 // give_nne ok 00044 // gt->give_nodes ok 00045 // gt->give_node_coord2d ok 00046 // */ 00047 // /** 00048 // */ 00049 // void newmeshread (const char *filename,long lcid) 00050 // { 00051 // long ip_dir,nod_spr; 00052 // ip_dir = 0; 00053 // nod_spr = 1; 00054 // 00055 // if ((Mp->adaptflag & 16)) 00056 // nod_spr = 0; 00057 // 00058 // 00059 // // ********************************************************************** 00060 // // begin of step 1 - aproximation `values` into nodes 00061 // // ********************************************************************** 00062 // 00063 // long ndofn,dim,nipcomp[3]; 00064 // double *nodvals_ip; 00065 // 00066 // dim = 2; 00067 // ndofn = 2; 00068 // 00069 // nipcomp[0] = 0; // number of transferd components of array ip[]->strain 00070 // nipcomp[1] = 0; // number of transferd components of array ip[]->stress 00071 // nipcomp[2] = 5; // number of transferd components of array ip[]->eqother !!! bez hardeningu 00072 // 00073 // if (!ip_dir) 00074 // give_nodvals_ip (dim,nipcomp,nodvals_ip); 00075 // else 00076 // nodvals_ip = NULL; 00077 // 00078 // // ********************************************************************** 00079 // // begin of step 2 - loading of new mesh 00080 // // ********************************************************************** 00081 // 00082 // problem *p_old = new problem; 00083 // problem *p_new = new problem; 00084 // 00085 // p_old->globinic (); 00086 // 00087 // Mp = NULL; Gtm = NULL; Mt = NULL; Mm = NULL; Mc = NULL; Mb = NULL; Lsrs = NULL; Smat = NULL; 00088 // fclose (Out); 00089 // 00090 // stochdriver stochd; 00091 // const char *argv[2]; 00092 // argv[0] = ""; 00093 // argv[1] = filename; 00094 // mefel_init (2, argv, &stochd); 00095 // 00096 // p_new->globinic (); 00097 // 00098 // // ********************************************************************** 00099 // // begin of step 3 - konverse 00100 // // ********************************************************************** 00101 // 00102 // long *parentel_ip,*parentel_nod; 00103 // double **ipcoord_new; 00104 // 00105 // parentel_ip = new long [p_new->mm->tnip]; 00106 // parentel_nod = new long [p_new->mt->nn]; 00107 // allocm (p_new->mm->tnip,3,ipcoord_new); 00108 // 00109 // 00110 // findout_parentel_ip (p_old->gt,Gtm,Mm,Mt,ipcoord_new,parentel_ip); 00111 // findout_parentel_nod (p_old->gt,Gtm,parentel_nod,dim); 00112 // 00113 // //if (Mespr==1) fprintf (stdout,"\n *** Transformation of meshes ***\n"); 00114 // 00115 // if (!ip_dir) transfvalues_ip_indirect (p_old->gt,Mm, ipcoord_new,parentel_ip,nipcomp,nodvals_ip ); 00116 // else transfvalues_ip_direct (p_old->gt,Mm,(const double **)ipcoord_new,parentel_ip,nipcomp,p_old->mm,dim); 00117 // 00118 // if (!nod_spr) transfvalues_nod (p_old,p_new,lcid,dim,ndofn,parentel_nod,'n'); 00119 // else transfvalues_nod (p_old,p_new,lcid,dim,ndofn,parentel_nod,'y'); 00120 // 00121 // 00122 // delete [] parentel_ip; 00123 // delete [] parentel_nod; 00124 // destrm (ipcoord_new,Mm->tnip); 00125 // 00126 // // ********************************************************************** 00127 // // end of step 3 - konverse 00128 // // ********************************************************************** 00129 // 00130 // p_old->dealoc (); 00131 // p_new->deinic (); 00132 // 00133 // // ********************************************************************** 00134 // // end of step 2 - loading of new mesh 00135 // // ********************************************************************** 00136 // 00137 // if (nodvals_ip!=NULL) delete [] nodvals_ip; 00138 // 00139 // // ********************************************************************** 00140 // // end of step 1 - aproximation `values` into nodes 00141 // // ********************************************************************** 00142 // 00143 // 00144 // mefel_right_hand_side (lcid,Lsrs->rhs); 00145 // } 00146 // 00147 // /** 00148 // Function approximates values stored on integration points (in arrays `strain` `stress` `eqother`) to nodes - for "global" mesh. 00149 // Used interpolation method is 'spr_smoothing'. 00150 // Values in nodes are returned by array nodvalues. 00151 // Value position in array is (val[0] on nod[0] ... val[0] on nod[gt->nn] ... val[nval] on nod[0] ... val[nval] on nod[gt->nn]) 00152 // Necessary precondition is Mm->ip[i].ncompstr and Mm->ip[i].ncompother are constant all over the domain !!!! 00153 // 00154 // @param dim - dimension of solved problem 00155 // @param nipcomp - number of members (of arrays Mm->ip[].strain , Mm->ip[].stress and Mm->ip[].eqother), which are approximated; 00156 // nipcomp[0] <0;Mm->ip[].ncompstr>, nipcomp[1] <0;Mm->ip[].ncompstr>, nipcomp[2] <0;Mm->ip[].ncompother> 00157 // if nipcomp = -1 => all members 00158 // @param nodvalues - answer ; nonalocated 00159 // 00160 // created 9.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00161 // */ 00162 // long give_nodvals_ip (long dim,long nipcomp[],double *&nodvalues) 00163 // { 00164 // long i,j,k,nip,ipp,nvals,*nsp; 00165 // double **spcoord,**spvalue,**ipcoord; 00166 // 00167 // if (nipcomp[0]==-1 || nipcomp[0] > Mm->ip[0].ncompstr ) nipcomp[0] = Mm->ip[0].ncompstr; 00168 // if (nipcomp[1]==-1 || nipcomp[1] > Mm->ip[0].ncompstr ) nipcomp[1] = Mm->ip[0].ncompstr; 00169 // if (nipcomp[2]==-1 || nipcomp[2] > Mm->ip[0].ncompother) nipcomp[2] = Mm->ip[0].ncompother; 00170 // 00171 // nvals = nipcomp[0] + nipcomp[1] + nipcomp[2]; 00172 // 00173 // if (nvals==0){ 00174 // nodvalues = NULL; 00175 // return 0; 00176 // } 00177 // else 00178 // nodvalues = new double [nvals*Mt->nn]; 00179 // 00180 // nsp = new long [Mt->ne]; 00181 // spcoord = new double* [Mt->ne]; 00182 // spvalue = new double* [Mt->ne]; 00183 // ipcoord = new double* [Mm->tnip]; 00184 // for (i=0;i<Mm->tnip;i++) 00185 // ipcoord[i] = new double [3]; 00186 // 00187 // allipcoord (ipcoord); 00188 // 00189 // for (i=0;i<Mt->ne;i++){ 00190 // nip = 0; 00191 // for (j=0;j<Mt->elements[i].nb;j++) 00192 // for (k=0;k<Mt->elements[i].nb;k++) 00193 // nip += Mt->elements[i].nip[j][k]; 00194 // 00195 // nsp[i] = nip; 00196 // 00197 // ipp = Mt->elements[i].ipp[0][0]; 00198 // spcoord[i] = new double [nip*dim]; 00199 // spvalue[i] = new double [nip*nvals]; 00200 // for (j=0;j<nip;j++){ 00201 // for (k=0;k<dim;k++) 00202 // spcoord[i][j*dim+k] = ipcoord[ipp+j][k]; 00203 // 00204 // for (k=0;k<nipcomp[0];k++) 00205 // spvalue[i][j*nvals + k ] = Mm->ip[ipp+j].strain[k]; 00206 // 00207 // for (k=0;k<nipcomp[1];k++) 00208 // spvalue[i][j*nvals + k + nipcomp[0] ] = Mm->ip[ipp+j].stress[k]; 00209 // 00210 // for (k=0;k<nipcomp[2];k++) 00211 // spvalue[i][j*nvals + k + nipcomp[0] + nipcomp[1]] = Mm->ip[ipp+j].eqother[k]; 00212 // 00213 // } 00214 // } 00215 // 00216 // least_square spr(dim,nvals,Gtm,0+((Mp->adaptflag & 16) ? 16 : 0)); // co to je 1 ?? 00217 // 00218 // spr.spr_optional (Out,nsp,spcoord,spvalue,nodvalues,1); 00219 // 00220 // delete [] nsp; 00221 // for (i=0;i<Mt->ne;i++){ delete [] spcoord[i]; delete [] spvalue[i]; } 00222 // delete [] spcoord; 00223 // delete [] spvalue; 00224 // for (i=0;i<Mm->tnip;i++) delete [] ipcoord[i]; 00225 // delete [] ipcoord; 00226 // 00227 // return nvals; 00228 // } 00229 // 00230 // /** 00231 // Function returns array of deformations in ALL nodes of problem 'p'. 00232 // Deformation position in answer array is (r[node_1][x] ... r[node_nn][x] , r[node_1][y] ... r[node_nn][y]) 00233 // 00234 // @param lcid - id of load case 00235 // @param p - pointr on problem 00236 // @param rfull - answer - allocated 1D array, size is p->mt->nn * p->mm->ip[0].ndofn 00237 // 00238 // created 9.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00239 // */ 00240 // void give_rfull (long lcid,problem *p,double *rfull) 00241 // { 00242 // long i,j,ndofn; 00243 // double *r; 00244 // 00245 // ndofn = p->mt->give_ndofn (0); 00246 // r = new double [ndofn]; 00247 // 00248 // for (i=0;i<p->mt->nn;i++){ 00249 // noddispl (lcid,i,p,r); 00250 // for (j=0;j<ndofn;j++) 00251 // rfull[i+j*p->mt->nn] = r[j]; 00252 // } 00253 // 00254 // delete [] r; 00255 // } 00256 // 00257 // /** 00258 // function computes coordinates of all integration points in "global" mesh 00259 // 00260 // @param ipcoord - empty 1D array, size is Mt->tnip x 3 00261 // 00262 // created 5.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00263 // */ 00264 // void allipcoord (double **ipcoord) 00265 // { 00266 // long i,ri,ci,nb,ipp; 00267 // 00268 // for (i=0;i<Mt->ne;i++){ 00269 // nb = Mt->elements[i].nb; 00270 // for (ri=0;ri<nb;ri++){ 00271 // for (ci=0;ci<nb;ci++){ 00272 // ipp = Mt->elements[i].ipp[ri][ci]; 00273 // 00274 // if (Mt->elements[i].nip[ri][ci]){ 00275 // switch (Mt->give_elem_type(i)){ 00276 // case planeelementlt:{ Pelt->ipcoordblock (i,ri,ci,ipcoord+ipp); break; } 00277 // case planeelementqt:{ Peqt->ipcoordblock (i,ri,ci,ipcoord+ipp); break; } 00278 // case planeelementlq:{ Pelq->ipcoordblock (i,ri,ci,ipcoord+ipp); break; } 00279 // case planeelementqq:{ Peqq->ipcoordblock (i,ri,ci,ipcoord+ipp); break; } 00280 // default:{ 00281 // fprintf (stderr,"\n\n unknown element type is required in function allipcoord (%s, line %d).\n",__FILE__,__LINE__); 00282 // } 00283 // } 00284 // } 00285 // } 00286 // } 00287 // } 00288 // } 00289 // 00290 // /** 00291 // Function alocates required arrays for function 'locintpoints', 00292 // which finds out (for every integration point from new mesh) "parent" element (from old mesh), 00293 // in which the integration point lays. 00294 // Parent elements are returned in array "parentel_ip". 00295 // 00296 // @param gt_old - gtopology of old mesh 00297 // @param gt_new - gtopology of new mesh 00298 // @param mm_new - mechmat of new mesh 00299 // @param mt_new - mechtop of new mesh 00300 // @param ipcoord - array of coordinates of new mesh int. points, size is mm_new->tnip x dim 00301 // @param parentel_ip - answer = 1D array, size is mt_new->tnip 00302 // 00303 // created 3.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00304 // */ 00305 // void findout_parentel_ip (gtopology *gt_old,gtopology *gt_new,mechmat *mm_new,mechtop *mt_new,double **ipcoord,long *parentel_ip) 00306 // { 00307 // long *passedel1; passedel1 = new long [gt_old->ne]; 00308 // long *passedel2; passedel2 = new long [gt_new->ne]; memset (passedel2,0,gt_new->ne*sizeof(long)); 00309 // long *elheap; elheap = new long [gt_new->ne]; elheap[0] = 0; 00310 // long *newelheap; newelheap = new long [gt_new->ne]; 00311 // long *susel; susel = new long [gt_old->ne]; 00312 // long *newsusel; newsusel = new long [gt_old->ne]; 00313 // 00314 // memset (parentel_ip,0,mm_new->tnip*sizeof(long)); 00315 // allipcoord (ipcoord); 00316 // 00317 // locintpoints (gt_old,gt_new,mt_new,ipcoord,passedel1,passedel2,1,elheap,newelheap,susel,newsusel,parentel_ip); 00318 // 00319 // delete [] passedel1; 00320 // delete [] passedel2; 00321 // delete [] elheap; 00322 // delete [] newelheap; 00323 // delete [] susel; 00324 // delete [] newsusel; 00325 // } 00326 // 00327 // /** 00328 // For every integration point (from new mesh) function finds out "parent" element (from old mesh), 00329 // in which the int. point lays. 00330 // Parent elements are returned in array "parentel". 00331 // 00332 // @param gt_old - gtopology of old mesh 00333 // @param gt_new - gtopology of new mesh 00334 // @param mt_new - mechtop of new mesh 00335 // @param ipcoord - array of coordinates of new mesh int. points, size is mm_new->tnip x dim 00336 // @param passedel1 - empty array, size is gt_old->ne 00337 // @param passedel2 - zero array, size is gt_new->ne 00338 // @param nelheap = 1 00339 // @param elheap - array of size gt_new->ne, elheap[0] = 0 00340 // @param newelheap - empty array, size is gt_new->ne 00341 // @param susel - empty array, size is gt_old->ne 00342 // @param newsusel - answer = 1D array, size is gt_old->ne 00343 // @param parentel - answer = 1D array, size is mm_new->tnip 00344 // 00345 // created 3.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00346 // */ 00347 // void locintpoints (gtopology *gt_old,gtopology *gt_new,mechtop *mt_new,double **ipcoord,long *passedel1,long *passedel2,long nelheap,long *elheap,long *newelheap,long *susel,long *newsusel,long *parentel) 00348 // { 00349 // long i,j,k,l,el,nip,ipp,nnewelheap,nsusel; 00350 // 00351 // nnewelheap = 0; 00352 // for (i=0;i<nelheap;i++){ 00353 // nsusel = mt_new->give_nip(elheap[i],0,0); 00354 // ipp = mt_new->elements[elheap[i]].ipp[0][0]; 00355 // for (j=0;j<nsusel;j++) 00356 // susel[j] = parentel[ipp++]; 00357 // 00358 // for (j=0;j<gt_new->nadjelnod[elheap[i]];j++){ 00359 // el = gt_new->adjelel[elheap[i]][j]; 00360 // 00361 // if (!passedel2[el]){ 00362 // nip = 0; 00363 // ipp = mt_new->elements[el].nb; 00364 // for (k=0;k<ipp;k++) 00365 // for (l=0;l<ipp;l++) 00366 // nip += mt_new->elements[el].nip[k][l]; 00367 // 00368 // ipp = mt_new->elements[el].ipp[0][0]; 00369 // for (k=0;k<nip;k++){ 00370 // for (l=0;l<gt_old->ne;l++) 00371 // passedel1[l] = 0; 00372 // 00373 // parentel[ipp+k] = whereispoint (gt_old,nsusel,susel,newsusel,passedel1,ipcoord[ipp+k][0],ipcoord[ipp+k][1],'f'); 00374 // } 00375 // 00376 // passedel2[el] = 1; 00377 // newelheap[nnewelheap++] = el; 00378 // } 00379 // } 00380 // } 00381 // 00382 // if (!nnewelheap) 00383 // return ; 00384 // 00385 // locintpoints (gt_old,gt_new,mt_new,ipcoord,passedel1,passedel2,nnewelheap,newelheap,elheap,susel,newsusel,parentel); 00386 // return ; 00387 // } 00388 // 00389 // /** 00390 // Function runs function 'locnodes'. 00391 // 00392 // @param gt_old - gtopology of old mesh 00393 // @param gt_new - gtopology of new mesh 00394 // @param parentel_nod - answer = array of parent elements, size is mt_new->nn 00395 // 00396 // created 6.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00397 // */ 00398 // void findout_parentel_nod (gtopology *gt_old,gtopology *gt_new,long *parentel_nod,long dim) 00399 // { 00400 // locnodes (gt_old,gt_new,parentel_nod,dim); 00401 // } 00402 // 00403 // /** 00404 // For every node (from new mesh) function finds out "parent" element (from old mesh), 00405 // in which the "new" node lays. 00406 // 00407 // @param gt_old - gtopology of old mesh 00408 // @param gt_new - gtopology of new mesh 00409 // @param parentel - answer = allocated array of parent elements, size is mt_new->nn 00410 // 00411 // created 4.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00412 // */ 00413 // void locnodes (gtopology *gt_old,gtopology *gt_new,long *parentel,long dim) 00414 // { 00415 // long i,j,k,n,nod,endheap; 00416 // ivector oldsusel(gt_old->ne),newsusel(gt_old->ne),passedel(gt_old->ne),nodheap(gt_new->nn),nsusel(gt_new->nn); 00417 // imatrix susel; 00418 // 00419 // if (dim==2) allocm (gt_new->nn,18,susel); 00420 // if (dim==3) allocm (gt_new->nn,120,susel); 00421 // 00422 // fillv (-1,parentel,gt_new->nn); 00423 // 00424 // endheap = 1; 00425 // nsusel[0] = 1; 00426 // susel[0][0] = 0; 00427 // nodheap[0] = 0; 00428 // 00429 // for (i=0;i<gt_new->nn;i++){ 00430 // nod = nodheap[i]; 00431 // 00432 // nullv (passedel); 00433 // copyv (susel[nod],oldsusel.a,nsusel[nod]); 00434 // 00435 // parentel[nod] = whereispoint (gt_old,nsusel[nod],oldsusel.a,newsusel.a,passedel.a,gt_new->gnodes[nod].x,gt_new->gnodes[nod].y,'f'); 00436 // 00437 // for (j=0;j<gt_new->nadjelnod[nod];j++) 00438 // for (k=0;k<gt_new->gelements[gt_new->adjelnod[nod][j]].nne;k++){ 00439 // n = gt_new->gelements[gt_new->adjelnod[nod][j]].nodes[k]; 00440 // if (parentel[n] == -1){ 00441 // parentel[n] = -2; 00442 // nodheap[endheap++] = n; 00443 // } 00444 // if (parentel[n] == -2) 00445 // susel[n][nsusel[n]++] = parentel[nod]; 00446 // } 00447 // 00448 // } 00449 // } 00450 // 00451 // /** 00452 // Function finds out element including "point". 00453 // First of all, the elements in "susel" are investigated. 00454 // 00455 // @param gt - gtopology, where the element is finding out 00456 // @param nsusel - number suspected elements 00457 // @param susel - array of suspected elements, size is gt->ne 00458 // @param newsusel - empty array, size is gt->ne 00459 // @param passedel - zero !!!!! array, size is gt->ne 00460 // @param x,y - coordinates of point 00461 // 00462 // created 3.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00463 // */ 00464 // long whereispoint (gtopology *gt,long nsusel,long *susel,long *newsusel,long *passedel,double x,double y,char flag) 00465 // { 00466 // long i,j,el,nnewsusel; 00467 // 00468 // if (flag=='f') 00469 // for (i=0;i<nsusel;i++) 00470 // if ( !passedel[susel[i]]++ && (ispointinel (gt->gnodes,gt->gelements[susel[i]].nodes,x,y,gt->gelements[susel[i]].auxinf)) ) 00471 // return (susel[i]); 00472 // 00473 // nnewsusel = 0; 00474 // for (i=0;i<nsusel;i++){ 00475 // for (j=0;j<gt->nadjelnod[susel[i]];j++){ 00476 // el = gt->adjelel[susel[i]][j]; 00477 // if (!passedel[el]){ 00478 // if (ispointinel (gt->gnodes,gt->gelements[el].nodes,x,y,gt->gelements[el].auxinf)) 00479 // return (el); 00480 // 00481 // passedel[el] = 1; 00482 // newsusel[nnewsusel++] = el; 00483 // } 00484 // } 00485 // } 00486 // 00487 // if (!nnewsusel){ 00488 // fprintf (stderr,"\n\n!!!!!!!!! Point does not lay in domain (%s, line %d)\n\n",__FILE__,__LINE__); 00489 // return (-1); 00490 // } 00491 // 00492 // return ( whereispoint (gt,nnewsusel,newsusel,susel,passedel,x,y,'n') ); 00493 // } 00494 // 00495 // /** 00496 // Function finds out, whether 'node' lays in 'element'. 00497 // Node is defined by its coordinates - x,y. 00498 // Element is defined by its node numbers - 'elnod' and dimdegnne - 'ndd' 00499 // 00500 // @param nod - array of all gnodes 00501 // @param elnod - array of node numbers of element 00502 // @param x,y - coordinates of node 00503 // @param ndd - nnedimdeg of element 00504 // 00505 // created 3.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00506 // */ 00507 // long ispointinel (gnode *nod,long *elnod,double x,double y,long ndd) 00508 // { 00509 // switch (ndd){ 00510 // case 312:{ 00511 // if ( isnodonlhsofline (nod[elnod[0]],nod[elnod[1]],x,y) && 00512 // isnodonlhsofline (nod[elnod[1]],nod[elnod[2]],x,y) && 00513 // isnodonlhsofline (nod[elnod[2]],nod[elnod[0]],x,y) ) 00514 // return (1); 00515 // break; 00516 // } 00517 // case 412:{ 00518 // if ( isnodonlhsofline (nod[elnod[0]],nod[elnod[1]],x,y) && 00519 // isnodonlhsofline (nod[elnod[1]],nod[elnod[2]],x,y) && 00520 // isnodonlhsofline (nod[elnod[2]],nod[elnod[3]],x,y) && 00521 // isnodonlhsofline (nod[elnod[3]],nod[elnod[0]],x,y) ) 00522 // return (1); 00523 // break; 00524 // } 00525 // case 622:{ 00526 // if ( isnodonlhsof3pcurve (nod[elnod[0]],nod[elnod[3]],nod[elnod[1]],x,y) && 00527 // isnodonlhsof3pcurve (nod[elnod[1]],nod[elnod[4]],nod[elnod[2]],x,y) && 00528 // isnodonlhsof3pcurve (nod[elnod[2]],nod[elnod[5]],nod[elnod[0]],x,y) ) 00529 // return (1); 00530 // break; 00531 // } 00532 // case 822:{ 00533 // if ( isnodonlhsof3pcurve (nod[elnod[0]],nod[elnod[4]],nod[elnod[1]],x,y) && 00534 // isnodonlhsof3pcurve (nod[elnod[1]],nod[elnod[5]],nod[elnod[2]],x,y) && 00535 // isnodonlhsof3pcurve (nod[elnod[2]],nod[elnod[6]],nod[elnod[3]],x,y) && 00536 // isnodonlhsof3pcurve (nod[elnod[3]],nod[elnod[7]],nod[elnod[0]],x,y) ) 00537 // return (1); 00538 // break; 00539 // } 00540 // //case 413:{ break; } 00541 // default:{ 00542 // fprintf (stderr,"\n\n wrong dimdegnne in function ispointinel (%s, line %d)",__FILE__,__LINE__); 00543 // break; 00544 // }} 00545 // 00546 // return (0); 00547 // } 00548 // 00549 // /** 00550 // Function finds out, whether point (xx,yy) lays on the left hand side of 'line'. 00551 // Line is defined by two nodes n1,n2. 00552 // 00553 // created 3.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00554 // */ 00555 // long isnodonlhsofline (gnode &n1,gnode &n2,double xx,double yy) 00556 // { 00557 // if (n1.x != n2.x) 00558 // if ( (n2.x>n1.x ? 1.0:-1.0)*((xx-n1.x)*(n1.y-n2.y)/(n1.x-n2.x)+(n1.y-yy)) <= 1.0e-13 ) 00559 // return (1); 00560 // else 00561 // return (0); 00562 // else 00563 // if ( (n1.y-n2.y)*(xx-n1.x) > -1.0e-13 ) 00564 // return (1); 00565 // else 00566 // return (0); 00567 // } 00568 // 00569 // /** 00570 // Function finds out, whether point (xx,yy) lays on the left hand side of 'curve'. 00571 // Curve is defined by three nodes n1,n2,n3. 00572 // 00573 // created 3.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00574 // */ 00575 // long isnodonlhsof3pcurve (gnode &n1,gnode &n2,gnode &n3,double xx,double yy) 00576 // { 00577 // double c; 00578 // 00579 // c = (n1.x*(n3.y-n2.y)+n2.x*(n1.y-n3.y)+n3.x*(n2.y-n1.y))/(fabs(n1.x-n3.x)+fabs(n1.y-n3.y)); 00580 // if (-1.0e-13<c && c<1.0e-13) 00581 // return (isnodonlhsofline (n1,n3,xx,yy)); 00582 // else{ 00583 // // docasne 00584 // //fprintf (stdout,"\n\n\n hrana je krivocara \n\n\n"); 00585 // // docasne 00586 // double s,x[3],y[2]; 00587 // 00588 // c = sqrt((n3.x-n1.x)*(n3.x-n1.x)+(n3.y-n1.y)*(n3.y-n1.y)); 00589 // s = (n3.y-n1.y)/c; 00590 // c = (n3.x-n1.x)/c; 00591 // 00592 // x[0] = c*( xx-n1.x) + s*( yy-n1.y); y[0] = -s*( xx-n1.x) + c*( yy-n1.y); 00593 // x[1] = c*(n2.x-n1.x) + s*(n2.y-n1.y); y[1] = -s*(n2.x-n1.x) + c*(n2.y-n1.y); 00594 // x[2] = c*(n3.x-n1.x) + s*(n3.y-n1.y); 00595 // 00596 // if ( y[1]*x[0]*x[0]/x[1]/(x[1]-x[2]) + y[1]*x[2]*x[0]/x[1]/(x[2]-x[1]) -y[0] <= 1.0e-13 ) 00597 // return (1); 00598 // 00599 // return (0); 00600 // } 00601 // } 00602 // 00603 // /** 00604 // Function transforms values from nodes of old mesh into int. points of new mesh. 00605 // 00606 // @param gt_old - gtopology of old mesh 00607 // @param mm_new - mechmat of new mesh 00608 // @param ipcoord - array of coordinates of new mesh int. points, size is mm_new->tnip x dim 00609 // @param parentel - array of parent elements of new mesh int. points, size is mt_new->tnip 00610 // @param nipcomp - number of members (of arrays Mm->ip[].strain , Mm->ip[].stress and Mm->ip[].eqother), which are approximated; 00611 // nipcomp[0] <0;Mm->ip[].ncompstr>, nipcomp[1] <0;Mm->ip[].ncompstr>, nipcomp[2] <0;Mm->ip[].ncompother> 00612 // if nipcomp = -1 => all members 00613 // @param nodvalue - 1D array of nodal values of old mesh, size is gt_old->nn * nvals 00614 // value position in array is (val[0] on nod[0] ... val[0] on nod[gt->nn] ... val[nval] on nod[0] ... val[nval] on nod[gt->nn]) 00615 // 00616 // created 7.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00617 // */ 00618 // void transfvalues_ip_indirect (gtopology *gt_old,mechmat *mm_new,double **ipcoord,long *parentel,long nipcomp[],double *nodvalue) 00619 // { 00620 // long i,j,maxnchild,nvals; 00621 // long *nchilds,**childip; 00622 // double *x,*y,**pointval1,**pointval2,**pointval3; 00623 // 00624 // 00625 // if (nipcomp[0]==-1 || nipcomp[0] > Mm->ip[0].ncompstr ) nipcomp[0] = Mm->ip[0].ncompstr; 00626 // if (nipcomp[1]==-1 || nipcomp[1] > Mm->ip[0].ncompstr ) nipcomp[1] = Mm->ip[0].ncompstr; 00627 // if (nipcomp[2]==-1 || nipcomp[2] > Mm->ip[0].ncompother) nipcomp[2] = Mm->ip[0].ncompother; 00628 // 00629 // nvals = nipcomp[0] + nipcomp[1] + nipcomp[2]; 00630 // 00631 // nchilds = new long [gt_old->ne]; 00632 // memset (nchilds,0,gt_old->ne*sizeof(long)); 00633 // 00634 // for (i=0;i<mm_new->tnip;i++) 00635 // nchilds[parentel[i]]++; 00636 // 00637 // maxnchild = 0; 00638 // childip = new long* [gt_old->ne]; 00639 // for (i=0;i<gt_old->ne;i++){ 00640 // childip[i] = new long [nchilds[i]]; 00641 // if (maxnchild<nchilds[i]) maxnchild = nchilds[i]; 00642 // nchilds[i] = 0; 00643 // } 00644 // 00645 // for (i=0;i<mm_new->tnip;i++) 00646 // childip[parentel[i]][nchilds[parentel[i]]++] = i; 00647 // 00648 // x = new double [maxnchild]; 00649 // y = new double [maxnchild]; 00650 // pointval1 = new double* [maxnchild]; 00651 // pointval2 = new double* [maxnchild]; 00652 // pointval3 = new double* [maxnchild]; 00653 // 00654 // for (i=0;i<gt_old->ne;i++){ 00655 // for (j=0;j<nchilds[i];j++){ 00656 // x[j] = ipcoord[childip[i][j]][0]; 00657 // y[j] = ipcoord[childip[i][j]][1]; 00658 // 00659 // if (nipcomp[0]) pointval1[j] = mm_new->ip[childip[i][j]].strain; 00660 // if (nipcomp[1]) pointval2[j] = mm_new->ip[childip[i][j]].stress; 00661 // if (nipcomp[2]) pointval3[j] = mm_new->ip[childip[i][j]].eqother; 00662 // } 00663 // 00664 // if (nipcomp[0]) give_valuesinpoints (gt_old,i,nchilds[i],x,y,nipcomp[0],nodvalue ,pointval1,'2'); 00665 // if (nipcomp[1]) give_valuesinpoints (gt_old,i,nchilds[i],x,y,nipcomp[1],nodvalue + gt_old->nn*(nipcomp[0] ),pointval2,'2'); 00666 // if (nipcomp[2]) give_valuesinpoints (gt_old,i,nchilds[i],x,y,nipcomp[2],nodvalue + gt_old->nn*(nipcomp[0]+nipcomp[1]),pointval3,'2'); 00667 // } 00668 // 00669 // 00670 // delete [] nchilds; 00671 // for (i=0;i<gt_old->ne;i++) delete [] childip[i]; 00672 // delete [] childip; 00673 // delete [] x; 00674 // delete [] y; 00675 // for (i=0;i<maxnchild;i++) pointval1[i] = NULL; delete [] pointval1; 00676 // for (i=0;i<maxnchild;i++) pointval2[i] = NULL; delete [] pointval2; 00677 // for (i=0;i<maxnchild;i++) pointval3[i] = NULL; delete [] pointval3; 00678 // } 00679 // 00680 // /** 00681 // Function transforms values from int. points of old mesh into int. points of new mesh. 00682 // 00683 // @param gt_old - gtopology of old mesh 00684 // @param mm_new - mechmat of new mesh 00685 // @param apcoord - array of coordinates of new mesh int. points, size is mm_new->tnip x dim 00686 // @param parentel - array of parent elements of new mesh int. points, size is mt_new->tnip 00687 // @param nipcomp - number of members (of arrays Mm->ip[].strain , Mm->ip[].stress and Mm->ip[].eqother), which are approximated; 00688 // nipcomp[0] <0;Mm->ip[].ncompstr>, nipcomp[1] <0;Mm->ip[].ncompstr>, nipcomp[2] <0;Mm->ip[].ncompother> 00689 // if nipcomp = -1 => all members 00690 // @param mm_old - mechmat of old mesh 00691 // @param dim - dimension of solved problem 00692 // 00693 // created 7.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00694 // */ 00695 // void transfvalues_ip_direct (gtopology *gt_old,mechmat *mm_new,const double **apcoord,const long *parentel,long nipcomp[],const mechmat *mm_old,long dim) 00696 // { 00697 // long i,j,nvals; 00698 // double **spvalue,**apvalue; 00699 // 00700 // if (nipcomp[0]==-1 || nipcomp[0] > Mm->ip[0].ncompstr ) nipcomp[0] = Mm->ip[0].ncompstr; 00701 // if (nipcomp[1]==-1 || nipcomp[1] > Mm->ip[0].ncompstr ) nipcomp[1] = Mm->ip[0].ncompstr; 00702 // if (nipcomp[2]==-1 || nipcomp[2] > Mm->ip[0].ncompother) nipcomp[2] = Mm->ip[0].ncompother; 00703 // 00704 // nvals = nipcomp[0] + nipcomp[1] + nipcomp[2]; 00705 // 00706 // spvalue = new double* [gt_old->ne]; 00707 // 00708 // for (i=0;i<gt_old->ne;i++){ 00709 // switch (gt_old->gelements[i].auxinf){ 00710 // case 312:{ 00711 // spvalue[i] = new double [1*nvals]; 00712 // 00713 // for (j=0;j<nipcomp[0];j++) spvalue[i][j ] = mm_old->ip[i].strain[j]; 00714 // for (j=0;j<nipcomp[1];j++) spvalue[i][j+nipcomp[0] ] = mm_old->ip[i].stress[j]; 00715 // for (j=0;j<nipcomp[2];j++) spvalue[i][j+nipcomp[0]+nipcomp[1]] = mm_old->ip[i].eqother[j]; 00716 // 00717 // break; 00718 // } 00719 // default:{ 00720 // fprintf (stderr,"\n\n wrong dimdegnne in function transfvalues_ip_direct (%s, line %d)",__FILE__,__LINE__); 00721 // break; 00722 // }} 00723 // } 00724 // 00725 // //jede jen pro 312 !!!!!!!!!!!!!!!!!!!! 00726 // apvalue = new double* [mm_new->tnip]; 00727 // for (i=0;i<mm_new->tnip;i++) 00728 // apvalue[i] = mm_new->ip[i].eqother; // !!! zde to nereflektuje styl nstra x nstre x nother 00729 // 00730 // least_square spr(dim,nvals,gt_old,((Mp->adaptflag & 16) ? 16 : 0)); 00731 // spr.L2_sp2sp (Out,spvalue,mm_new->tnip,parentel,apcoord,apvalue); 00732 // 00733 // 00734 // destrm (spvalue,gt_old->ne); 00735 // for (i=0;i<mm_new->tnip;i++) apvalue[i] = NULL; 00736 // delete [] apvalue; 00737 // } 00738 // 00739 // /** 00740 // Function transforms displacement from nodes of old mesh into nodes of new mesh by direct interpolation. 00741 // 00742 // @param po - problem of old mesh 00743 // @param pn - problem of new mesh 00744 // @param lcid - id of load case 00745 // @param ndofn - number of DOFs in node 00746 // @param parentel - array of parent (old)elements of new mesh nodes, size is p_new->mt->nn 00747 // 00748 // created 7.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00749 // */ 00750 // //void transfvalues_nod (problem *po,problem *pn,long lcid,long ndofn,long *parentel) 00751 // 00752 // 00753 // /** 00754 // Function transforms values from nodes of old mesh into nodes of new mesh. 00755 // 00756 // @param po - problem of old mesh 00757 // @param pn - problem of new mesh 00758 // @param lcid - id of load case 00759 // @param dim - dimension of solved problem 00760 // @param ndofn - number of DOFs in node 00761 // @param parentel - array of parent (old)elements of new mesh nodes, size is p_new->mt->nn 00762 // 00763 // created 7.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00764 // */ 00765 // void transfvalues_nod (problem *po,problem *pn,long lcid,long dim,long ndofn,long *parentel,char spr) 00766 // { 00767 // long i,j,k,nid,maxnchild; 00768 // long *nchilds,**childnod; 00769 // double **spcoord, *r_old, **r_new; 00770 // vector x,y; 00771 // 00772 // // array of deformations of old mesh 00773 // r_old = new double [ndofn * po->gt->nn]; 00774 // give_rfull (lcid,po,r_old); 00775 // 00776 // // allocation 00777 // nchilds = new long [po->gt->ne]; 00778 // memset (nchilds,0,po->gt->ne*sizeof(long)); 00779 // 00780 // for (i=0;i<pn->gt->nn;i++) 00781 // nchilds[parentel[i]]++; 00782 // 00783 // childnod = new long* [po->gt->ne]; 00784 // r_new = new double* [po->gt->ne]; 00785 // for (i=0;i<po->gt->ne;i++){ 00786 // childnod[i] = new long [nchilds[i]]; 00787 // r_new[i] = new double [nchilds[i]*ndofn]; 00788 // } 00789 // 00790 // if (spr=='y'){ 00791 // spcoord = new double* [po->gt->ne]; 00792 // for (i=0;i<po->gt->ne;i++) spcoord[i] = new double [nchilds[i]*dim]; 00793 // } 00794 // else if (spr=='n'){ 00795 // maxnchild = 0; 00796 // for (i=0;i<po->gt->ne;i++) if (maxnchild<nchilds[i]) maxnchild = nchilds[i]; 00797 // allocv (maxnchild,x); allocv (maxnchild,y); 00798 // } 00799 // else fprintf (stderr,"\n\n unknown spr flag is required in function transfvalues_nod (%s, line %d).\n",__FILE__,__LINE__); 00800 // 00801 // // transformation from r_old to r_new 00802 // memset (nchilds,0,po->gt->ne*sizeof(long)); 00803 // for (i=0;i<pn->gt->nn;i++) 00804 // childnod[parentel[i]][nchilds[parentel[i]]++] = i; 00805 // 00806 // // ******************** 00807 // /* 00808 // double X,Y,f; 00809 // 00810 // for (i=0;i<po->mt->nn;i++){ 00811 // X = po->gt->gnodes[i].x; 00812 // Y = po->gt->gnodes[i].y; 00813 // f = -X*X*X*X/1000.0 - X*X*X/10.0 + X*X + 10*X -Y*Y*Y*Y/3.0 - Y*Y*Y/2.0 + 20*Y*Y + 30*Y - 212; 00814 // 00815 // r_old[i ] = f; 00816 // r_old[i+po->mt->nn] = f; 00817 // } 00818 // 00819 // long hod,min; 00820 // double sec = clock(); 00821 // */ 00822 // // ******************** 00823 // 00824 // if (spr=='y'){ 00825 // for (i=0;i<po->gt->ne;i++){ 00826 // for (j=0;j<nchilds[i];j++){ 00827 // spcoord[i][j*dim ] = pn->gt->gnodes[childnod[i][j]].x; 00828 // spcoord[i][j*dim+1] = pn->gt->gnodes[childnod[i][j]].y; 00829 // if (dim==3) spcoord[i][j*dim+2] = pn->gt->gnodes[childnod[i][j]].z; 00830 // } 00831 // } 00832 // 00833 // least_square spr(dim,ndofn,po->gt,((Mp->adaptflag & 16) ? 16 : 0)); 00834 // spr.L2_nod2sp (Out,nchilds,spcoord,r_old,r_new,'n'); 00835 // } 00836 // else if (spr=='n'){ 00837 // for (i=0;i<po->gt->ne;i++){ 00838 // for (j=0;j<nchilds[i];j++){ 00839 // x[j] = pn->gt->gnodes[childnod[i][j]].x; 00840 // y[j] = pn->gt->gnodes[childnod[i][j]].y; 00841 // } 00842 // 00843 // give_valuesinpoints (po->gt,i,nchilds[i],x.a,y.a,ndofn,r_old,&r_new[i],'1'); 00844 // } 00845 // } 00846 // else fprintf (stderr,"\n\n unknown spr flag is required in function transfvalues_nod (%s, line %d).\n",__FILE__,__LINE__); 00847 // 00848 // // ******************** 00849 // /* 00850 // sec = (clock() - sec) / (double)CLOCKS_PER_SEC; 00851 // hod = (long)sec/3600; sec -= hod*3600; 00852 // min = (long)sec/60; sec -= min*60; 00853 // fprintf (stdout,"\n -----------------------------------"); 00854 // fprintf (stdout,"\n Consumed time by TRANSF %2ld:%02ld:%05.2f",hod,min,sec); 00855 // fprintf (stdout,"\n -----------------------------------\n"); 00856 // 00857 // double a,ai,nordr,norr=0.0; 00858 // vector dr(2*pn->gt->nn),r(pn->gt->nn); 00859 // 00860 // for (i=0;i<po->gt->ne;i++) 00861 // for (j=0;j<nchilds[i];j++){ 00862 // nid = childnod[i][j]; 00863 // X = pn->gt->gnodes[nid].x; 00864 // Y = pn->gt->gnodes[nid].y; 00865 // r[nid] = -X*X*X*X/1000.0 - X*X*X/10.0 + X*X + 10*X -Y*Y*Y*Y/3.0 - Y*Y*Y/2.0 + 20*Y*Y + 30*Y - 212; 00866 // 00867 // dr[nid+pn->gt->nn] = r[nid] - r_new[i][j*2+1]; 00868 // dr[nid ] = r[nid] - r_new[i][j*2 ]; 00869 // } 00870 // 00871 // sizev (dr,nordr); 00872 // fprintf (stdout,"\n\n nordr = %25.15f \n\n",nordr); 00873 // fprintf (stdout,"\n\n nordr/nn1 = %25.15f \n\n",nordr/sqrt(769)); 00874 // fprintf (stdout,"\n\n nordr/nn2 = %25.15f \n\n",nordr/sqrt(24048)); 00875 // 00876 // nordr = a = 0.0; 00877 // for (i=0;i<pn->gt->ne;i++){ 00878 // nordr += Pelt->error(i,dr,ai); 00879 // a += ai; 00880 // } 00881 // 00882 // nordr = sqrt(nordr/a); 00883 // fprintf (stdout,"\n\n nordr(A) = %25.15f \n\n",nordr); 00884 // 00885 // exit (1); 00886 // 00887 // 00888 // for (i=0;i<pn->gt->nn;i++) norr += r[i]; 00889 // sizev (r,norr); norr *= 2; 00890 // */ 00891 // // ******************** 00892 // 00893 // // saving deformations from r_old to p_new->lsrs->lhs 00894 // for (i=0;i<po->gt->ne;i++) 00895 // for (j=0;j<nchilds[i];j++){ 00896 // nid = childnod[i][j]; 00897 // for (k=0;k<ndofn;k++) 00898 // if (pn->gt->gnodes[nid].cn[k] > 0) 00899 // pn->lsrs->lhs[lcid*pn->lsrs->ndof + pn->gt->gnodes[nid].cn[k] - 1] = r_new[i][j*ndofn+k]; 00900 // } 00901 // 00902 // delete [] r_old; 00903 // delete [] nchilds; 00904 // destrm (childnod,po->gt->ne); 00905 // destrm (r_new,po->gt->ne); 00906 // if (spr=='y') destrm (spcoord,po->gt->ne); 00907 // } 00908 // 00909 // /** 00910 // Function approximates 'values' from nodes into 'points'. 00911 // All the point lays in 'element', which is determined by 'gt' and 'eid'. 00912 // Values in points are returned by array pointval. 00913 // 00914 // @param npoints - number of points 00915 // @param xx - array of x-coordinates of points, size is npoints 00916 // @param yy - array of y-coordinates of points, size is npoints 00917 // @param nval - number of values 00918 // @param nodvalues - 1D array of values in nodes, size is nval*gt->nn, value position in array is 00919 // (val[0] on nod[0] ... val[0] on nod[gt->nn] ... val[nval] on nod[0] ... val[nval] on nod[gt->nn]) 00920 // @param pointval - answer = 2D array, size is npoints*nval 00921 // 00922 // created 7.12.2002, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00923 // */ 00924 // void give_valuesinpoints (gtopology *gt,long eid,long npoints,double *xx,double *yy,long nval,double *nodvalues,double **pointval,char flag) 00925 // { 00926 // long i,j,k,nne; 00927 // double xi,eta,value; 00928 // vector x,y,nodval; 00929 // ivector nodes; 00930 // 00931 // nne = gt->give_nne (eid); 00932 // 00933 // allocv (nne,x); 00934 // allocv (nne,y); 00935 // allocv (nne,nodval); 00936 // allocv (nne,nodes); 00937 // 00938 // gt->give_nodes (eid,nodes); 00939 // gt->give_node_coord2d (x,y,eid); 00940 // 00941 // for (i=0;i<npoints;i++){ 00942 // 00943 // switch (gt->gelements[eid].auxinf){ 00944 // case 312:{ nc_lin_3_2d (xx[i],yy[i],x.a,y.a,xi,eta); break; } 00945 // case 622:{ nc_quad_3_2d (0.00000001,xx[i],yy[i],x.a,y.a,xi,eta); break; } 00946 // case 412:{ nc_lin_4_2d (0.00001,xx[i],yy[i],x.a,y.a,xi,eta); break; } 00947 // case 822:{ nc_quad_4_2d (0.00001,xx[i],yy[i],x.a,y.a,xi,eta); break; } 00948 // default:{ 00949 // fprintf (stderr,"\n\n unknown nnedegdim is required in function give_valuesinpoints (%s, line %d).\n",__FILE__,__LINE__); 00950 // }} 00951 // 00952 // for (j=0;j<nval;j++){ 00953 // for (k=0;k<nne;k++) 00954 // nodval[k] = nodvalues[nodes[k]+j*gt->nn]; 00955 // 00956 // switch (gt->gelements[eid].auxinf){ 00957 // case 312:{ value = Pelt->approx_nat (xi,eta,nodval); break; } 00958 // case 622:{ value = Peqt->approx (xi,eta,nodval); break; } 00959 // case 412:{ value = Pelq->approx (xi,eta,nodval); break; } 00960 // case 822:{ value = Peqq->approx (xi,eta,nodval); break; } 00961 // default:{ 00962 // fprintf (stderr,"\n\n unknown nnedegdim is required in function give_valuesinpoints (%s, line %d).\n",__FILE__,__LINE__); 00963 // }} 00964 // 00965 // if (flag=='2') 00966 // pointval[i][j] = value; 00967 // else if (flag=='1') 00968 // pointval[0][i*nval+j] = value; 00969 // else ; 00970 // } 00971 // } 00972 // } 00973 // 00974 /** 00975 function finds out node(from gt), which is the closest by required point 00976 00977 @param x,y,z - coordinates of point 00978 00979 created 30.1.2003, Ladislav Svoboda, termit@cml.fsv.cvut.cz 00980 */ 00981 long adjacnode (gtopology *gt,double x,double y,double z) 00982 { 00983 long i,a = 0; 00984 double d,min = 1e50; 00985 00986 for(i=0;i<gt->ne;i++){ 00987 d = (gt->gnodes[i].x - x)*(gt->gnodes[i].x - x) + (gt->gnodes[i].y - y)*(gt->gnodes[i].y - y) + (gt->gnodes[i].z - z)*(gt->gnodes[i].z - z); 00988 if (d<min){ 00989 min=d; 00990 a=i; 00991 } 00992 } 00993 return (a); 00994 } 00995 //////////////////// /* termitovo */ ////////////////////////////////////