# Contents of /branches/domexper/dudley/src/Util.c

Revision 765 - (show annotations)
Fri Jun 30 06:37:11 2006 UTC (13 years, 8 months ago) by gross
Original Path: trunk/finley/src/Util.c
File MIME type: text/plain
File size: 17732 byte(s)
The test for the contact normal has been modified to take in cosideration the fact that the normal is unique up to the factor +/-1.
Now the test checks the kllength of the normal for 1 and the angle to the reference normal.

 1 /* 2 ************************************************************ 3 * Copyright 2006 by ACcESS MNRF * 4 * * 5 * http://www.access.edu.au * 6 * Primary Business: Queensland, Australia * 7 * Licensed under the Open Software License version 3.0 * 8 * http://www.opensource.org/licenses/osl-3.0.php * 9 * * 10 ************************************************************ 11 */ 12 13 /**************************************************************/ 14 15 /* Some utility routines: */ 16 17 /**************************************************************/ 18 19 /* author: gross@access.edu.au */ 20 /* Version: \$Id\$ */ 21 22 /**************************************************************/ 23 24 #include "Finley.h" 25 #include "Util.h" 26 27 #ifdef _OPENMP 28 #include 29 #endif 30 31 /**************************************************************/ 32 33 /* returns true if any of the values in the short array values is not equalt to Zero */ 34 35 bool_t Finley_Util_anyNonZeroDouble(dim_t N, double* values) { 36 dim_t q; 37 for (q=0;q0) return TRUE; 38 return FALSE; 39 } 40 /**************************************************************/ 41 42 /* gathers double values out from in by index: */ 43 44 /* out(1:numData,1:len)=in(1:numData,index(1:len)) */ 45 46 void Finley_Util_Gather_double(dim_t len,index_t* index,dim_t numData,double* in, double * out){ 47 dim_t s,i; 48 for (s=0;s 0 ){ 151 det[q]=D; 152 D=1./D; 153 invA[q]=D; 154 } else { 155 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix"); 156 return; 157 } 158 } 159 break; 160 161 case 2: 162 for (q=0;q 0 ){ 170 det[q]=D; 171 D=1./D; 172 invA[INDEX3(0,0,q,2,2)]= A22*D; 173 invA[INDEX3(1,0,q,2,2)]=-A21*D; 174 invA[INDEX3(0,1,q,2,2)]=-A12*D; 175 invA[INDEX3(1,1,q,2,2)]= A11*D; 176 } else { 177 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix"); 178 return; 179 } 180 } 181 break; 182 183 case 3: 184 for (q=0;q 0 ){ 197 det[q] =D; 198 D=1./D; 199 invA[INDEX3(0,0,q,3,3)]=(A22*A33-A23*A32)*D; 200 invA[INDEX3(1,0,q,3,3)]=(A31*A23-A21*A33)*D; 201 invA[INDEX3(2,0,q,3,3)]=(A21*A32-A31*A22)*D; 202 invA[INDEX3(0,1,q,3,3)]=(A13*A32-A12*A33)*D; 203 invA[INDEX3(1,1,q,3,3)]=(A11*A33-A31*A13)*D; 204 invA[INDEX3(2,1,q,3,3)]=(A12*A31-A11*A32)*D; 205 invA[INDEX3(0,2,q,3,3)]=(A12*A23-A13*A22)*D; 206 invA[INDEX3(1,2,q,3,3)]=(A13*A21-A11*A23)*D; 207 invA[INDEX3(2,2,q,3,3)]=(A11*A22-A12*A21)*D; 208 } else { 209 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix"); 210 return; 211 } 212 } 213 break; 214 215 } 216 return; 217 } 218 219 /* sets the derterminate of a set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */ 220 221 void Finley_Util_DetOfSmallMat(dim_t len,dim_t dim,double* A, double* det){ 222 dim_t q; 223 register double A11,A12,A13,A21,A22,A23,A31,A32,A33; 224 225 switch(dim) { 226 case 1: 227 for (q=0;q0) { 279 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: area equals zero."); 280 return; 281 } else { 282 invlength=1./length; 283 Normal[INDEX2(0,q,2)]=A21*invlength; 284 Normal[INDEX2(1,q,2)]=-A11*invlength; 285 } 286 } 287 break; 288 case 3: 289 for (q=0;q0) { 301 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: area equals zero."); 302 return; 303 } else { 304 invlength=1./length; 305 Normal[INDEX2(0,q,3)]=CO_A13*invlength; 306 Normal[INDEX2(1,q,3)]=CO_A23*invlength; 307 Normal[INDEX2(2,q,3)]=CO_A33*invlength; 308 } 309 310 } 311 break; 312 313 } 314 return; 315 } 316 317 /* return the length of the vector which is orthogonal to the vectors A(:,0,q) and A(:,1,q) in the case of dim=3 */ 318 /* or the vector A(:,0,q) in the case of dim=2 */ 319 320 void Finley_LengthOfNormalVector(dim_t len, dim_t dim, dim_t dim1, double* A,double* length) { 321 dim_t q; 322 double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33; 323 324 switch(dim) { 325 case 1: 326 for (q=0;q=0) invMap[Map[i]]=i; 363 } 364 } 365 366 /* orders a Finley_Util_ValueAndIndex array by value */ 367 /* it is assumed that n is large */ 368 369 int Finley_Util_ValueAndIndex_compar(const void *arg1 , const void *arg2 ) { 370 Finley_Util_ValueAndIndex *e1,*e2; 371 e1=(Finley_Util_ValueAndIndex*) arg1; 372 e2=(Finley_Util_ValueAndIndex*) arg2; 373 if (e1->value < e2->value) return -1; 374 if (e1->value > e2->value) return 1; 375 if (e1->index < e2->index) return -1; 376 if (e1->index > e2->index) return 1; 377 378 return 0; 379 } 380 381 void Finley_Util_sortValueAndIndex(dim_t n,Finley_Util_ValueAndIndex* array) { 382 /* OMP : needs parallelization !*/ 383 qsort(array,n,sizeof(Finley_Util_ValueAndIndex),Finley_Util_ValueAndIndex_compar); 384 } 385 386 387 /**************************************************************/ 388 389 /* calculates the minimum value from a dim X N integer array */ 390 391 index_t Finley_Util_getMinInt(dim_t dim,dim_t N,index_t* values) { 392 dim_t i,j; 393 index_t out,out_local; 394 out=INDEX_T_MAX; 395 if (values!=NULL && dim*N>0 ) { 396 out=values[0]; 397 #pragma omp parallel private(out_local) 398 { 399 out_local=out; 400 #pragma omp for private(i,j) schedule(static) 401 for (j=0;j0 ) { 418 out=values[0]; 419 #pragma omp parallel private(out_local) 420 { 421 out_local=out; 422 #pragma omp for private(i,j) schedule(static) 423 for (j=0;j=0) { 441 index[out]=k; 442 out++; 443 } 444 } 445 return out; 446 } 447 448 /* returns true if array contains value */ 449 bool_t Finley_Util_isAny(dim_t N,index_t* array,index_t value) { 450 bool_t out=FALSE; 451 dim_t i; 452 #pragma omp parallel for private(i) schedule(static) reduction(||:out) 453 for (i=0;i=30 ) 514 fprintf( fid, "... " ); 515 fprintf( fid, "]\n" ); 516 } 517 void Finley_printIntArray( FILE *fid, dim_t n, int *array, char *name ) 518 { 519 index_t i; 520 521 if( name ) 522 fprintf( fid, "%s [ ", name ); 523 else 524 fprintf( fid, "[ " ); 525 for( i=0; i<(n<30 ? n : 30); i++ ) 526 fprintf( fid, "%d ", array[i] ); 527 if( n>=30 ) 528 fprintf( fid, "... " ); 529 fprintf( fid, "]\n" ); 530 } 531 void Finley_printMaskArray( FILE *fid, dim_t n, int *array, char *name ) 532 { 533 index_t i; 534 535 if( name ) 536 fprintf( fid, "%s [ ", name ); 537 else 538 fprintf( fid, "[ " ); 539 for( i=0; i<(n<30 ? n : 30); i++ ) 540 if( array[i]!=-1 ) 541 fprintf( fid, "%d ", array[i] ); 542 else 543 fprintf( fid, "* " ); 544 if( n>=30 ) 545 fprintf( fid, "... " ); 546 fprintf( fid, "]\n" ); 547 } 548 #endif 549 550 /* 551 * Revision 1.8 2005/08/12 01:45:43 jgs 552 * 553 * Revision 1.7.2.2 2005/09/07 06:26:22 gross 554 * the solver from finley are put into the standalone package paso now 555 * 556 * Revision 1.7.2.1 2005/08/04 22:41:11 gross 557 * some extra routines for finley that might speed-up RHS assembling in some cases (not actived right now) 558 * 559 * Revision 1.7 2005/07/08 04:07:59 jgs 560 * Merge of development branch back to main trunk on 2005-07-08 561 * 562 * Revision 1.1.1.1.2.4 2005/06/29 02:34:57 gross 563 * some changes towards 64 integers in finley 564 * 565 * Revision 1.1.1.1.2.3 2005/03/02 23:35:06 gross 566 * reimplementation of the ILU in Finley. block size>1 still needs some testing 567 * 568 * Revision 1.1.1.1.2.2 2005/02/18 02:27:31 gross 569 * two function that will be used for a reimplementation of the ILU preconditioner 570 * 571 * Revision 1.1.1.1.2.1 2004/11/12 06:58:19 gross 572 * a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry 573 * 574 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs 575 * initial import of project esys2 576 * 577 * Revision 1.3 2004/08/26 12:03:52 gross 578 * Some other bug in Finley_Assemble_gradient fixed. 579 * 580 * Revision 1.2 2004/07/02 04:21:13 gross 581 * Finley C code has been included 582 * 583 * Revision 1.1.1.1 2004/06/24 04:00:40 johng 584 * Initial version of eys using boost-python. 585 * 586 * 587 */

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision