Changeset 5821
- Timestamp:
- 09/15/10 10:18:38 (15 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Tria.cpp
r5803 r5821 1275 1275 1276 1276 int i; 1277 bool found = false;1278 1277 Input *input = NULL; 1279 1278 … … 2395 2394 void Tria::CreateKMatrixBalancedthickness_CG(Mat Kgg){ 2396 2395 2397 /* local declarations */ 2398 int i,j; 2399 2400 /* node data: */ 2396 /*Constants*/ 2401 2397 const int numdof=NDOF1*NUMVERTICES; 2402 double xyz_list[NUMVERTICES][3]; 2403 int* doflist=NULL; 2404 2405 /* gaussian points: */ 2406 int ig; 2407 GaussTria *gauss=NULL; 2408 2409 /* matrices: */ 2410 double L[NUMVERTICES]; 2411 double B[2][NUMVERTICES]; 2412 double Bprime[2][NUMVERTICES]; 2413 double DL[2][2]={0.0}; 2414 double DLprime[2][2]={0.0}; 2415 double DL_scalar; 2416 double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix 2417 double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. 2418 double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. 2419 double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. 2420 2421 double Jdettria; 2422 2423 /*input parameters for structural analysis (diagnostic): */ 2424 double dvx[2]; 2425 double dvy[2]; 2426 double vx,vy; 2427 double dvxdx,dvydy; 2428 double v_gauss[2]={0.0}; 2429 2430 double K[2][2]={0.0}; 2431 double KDL[2][2]={0.0}; 2432 int found=0; 2433 int dim; 2434 2435 /*inputs: */ 2436 Input* vxaverage_input=NULL; 2437 Input* vyaverage_input=NULL; 2438 bool artdiff; 2398 2399 /*Intermediaries */ 2400 bool artdiff; 2401 int i,j,ig,dim; 2402 double Jdettria ,vx,vy,dvxdx,dvydy; 2403 double dvx[2],dvy[2]; 2404 double xyz_list[NUMVERTICES][3]; 2405 double L[NUMVERTICES]; 2406 double B[2][NUMVERTICES]; 2407 double Bprime[2][NUMVERTICES]; 2408 double K[2][2] = {0.0}; 2409 double KDL[2][2] = {0.0}; 2410 double DL[2][2] = {0.0}; 2411 double DLprime[2][2] = {0.0}; 2412 double DL_scalar; 2413 double Ke_gg[numdof][numdof] = {0.0}; 2414 double Ke_gg_gaussian[numdof][numdof] = {0.0}; 2415 double Ke_gg_thickness1[numdof][numdof] = {0.0}; 2416 double Ke_gg_thickness2[numdof][numdof] = {0.0}; 2417 int *doflist = NULL; 2418 GaussTria *gauss = NULL; 2439 2419 2440 2420 /* Get node coordinates and dof list: */ … … 2442 2422 GetDofList(&doflist,NoneApproximationEnum,GsetEnum); 2443 2423 2444 /* retrieve someparameters: */2424 /*Retrieve all Inputs and parameters: */ 2445 2425 this->parameters->FindParam(&artdiff,ArtDiffEnum); 2446 2426 this->parameters->FindParam(&dim,DimEnum); 2447 2448 /*Retrieve all inputs we will be needed*/2427 Input* vxaverage_input=NULL; 2428 Input* vyaverage_input=NULL; 2449 2429 if(dim==2){ 2450 2430 vxaverage_input=inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input); … … 2458 2438 //Create Artificial diffusivity once for all if requested 2459 2439 if(artdiff){ 2460 //Get the Jacobian determinant2461 2440 gauss=new GaussTria(); 2462 2441 gauss->GaussCenter(); … … 2464 2443 delete gauss; 2465 2444 2466 //Build K matrix (artificial diffusivity matrix) 2467 vxaverage_input->GetParameterAverage(&v_gauss[0]); 2468 vyaverage_input->GetParameterAverage(&v_gauss[1]); 2469 2470 K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); 2471 K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); 2445 vxaverage_input->GetParameterAverage(&vx); 2446 vyaverage_input->GetParameterAverage(&vy); 2447 K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(vx); 2448 K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(vy); 2472 2449 } 2473 2450 … … 2478 2455 gauss->GaussPoint(ig); 2479 2456 2480 /* Get Jacobian determinant: */2481 2457 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); 2482 2483 /*Get B and B prime matrix: */2484 2458 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss); 2485 2459 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); 2486 2460 2487 //Get vx, vy and their derivatives at gauss point2488 2461 vxaverage_input->GetParameterValue(&vx,gauss); 2489 2462 vyaverage_input->GetParameterValue(&vy,gauss); 2490 2491 2463 vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); 2492 2464 vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); … … 2494 2466 dvxdx=dvx[0]; 2495 2467 dvydy=dvy[1]; 2496 2497 2468 DL_scalar=gauss->weight*Jdettria; 2498 2469 2499 //Create DL and DLprime matrix2500 2470 DL[0][0]=DL_scalar*dvxdx; 2501 2471 DL[1][1]=DL_scalar*dvydy; … … 2503 2473 DLprime[0][0]=DL_scalar*vx; 2504 2474 DLprime[1][1]=DL_scalar*vy; 2505 2506 //Do the triple product tL*D*L.2507 //Ke_gg_thickness=B'*DLprime*Bprime;2508 2475 2509 2476 TripleMultiply( &B[0][0],2,numdof,1, … … 2517 2484 &Ke_gg_thickness2[0][0],0); 2518 2485 2519 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */2520 2486 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness1[i][j]; 2521 2487 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_thickness2[i][j]; 2522 2488 2523 2489 if(artdiff){ 2524 2525 /* Compute artificial diffusivity */2526 2490 KDL[0][0]=DL_scalar*K[0][0]; 2527 2491 KDL[1][1]=DL_scalar*K[1][1]; … … 2532 2496 &Ke_gg_gaussian[0][0],0); 2533 2497 2534 /* Add artificial diffusivity matrix */2535 2498 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j]; 2536 2537 2499 } 2538 2500 } 2539 2501 2540 /*Add Ke_gg to global matrix Kgg: */2541 2502 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); 2542 2503 … … 2573 2534 /*input parameters for structural analysis (diagnostic): */ 2574 2535 double vx,vy; 2575 int found;2576 2536 int dim; 2577 2537 … … 2666 2626 double K[2][2]={0.0}; 2667 2627 double KDL[2][2]={0.0}; 2668 int found=0;2669 2628 int dim; 2670 2629 … … 2796 2755 void Tria::CreateKMatrixDiagnosticMacAyeal(Mat Kgg,Mat Kff, Mat Kfs){ 2797 2756 2798 const int numdof = 2 *NUMVERTICES;2757 /*Intermediaries*/ 2799 2758 double *Ke_gg_viscous = NULL; 2800 2759 double *Ke_gg_friction = NULL; … … 2805 2764 Ke_gg_friction = CreateKMatrixDiagnosticMacAyealFriction(); 2806 2765 2807 /* Initialize elementmatrix: */2766 /*Add Ke_gg values to Ke element stifness matrix: */ 2808 2767 Ke=this->NewElementMatrix(MacAyealApproximationEnum); 2809 2810 /*Add Ke_gg values to Ke element stifness matrix: */2811 2768 if(Ke_gg_viscous) Ke->AddValues(Ke_gg_viscous); 2812 2769 if(Ke_gg_friction)Ke->AddValues(Ke_gg_friction); … … 2819 2776 xfree((void**)&Ke_gg_viscous); 2820 2777 xfree((void**)&Ke_gg_friction); 2821 2822 } 2823 2778 } 2824 2779 /*}}}*/ 2825 2780 /*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealViscous{{{1*/ … … 3547 3502 double K[2][2]={0.0}; 3548 3503 double KDL[2][2]={0.0}; 3549 int found;3550 3504 int dim; 3551 3505 … … 3708 3662 /*input parameters for structural analysis (diagnostic): */ 3709 3663 double vx,vy; 3710 int found;3711 3664 int dim; 3712 3665 … … 3838 3791 3839 3792 int i,j; 3840 int found=0;3841 3793 3842 3794 /* node data: */ … … 5098 5050 void Tria::CreatePVectorThermalShelf( Vec pg){ 5099 5051 5100 int i ,found;5052 int i; 5101 5053 5102 5054 const int numdof=NUMVERTICES*NDOF1; … … 5186 5138 void Tria::CreatePVectorThermalSheet( Vec pg){ 5187 5139 5188 int i ,found;5140 int i; 5189 5141 5190 5142 const int numdof=NUMVERTICES*NDOF1; … … 6016 5968 int* sinternaldoflist=NULL; 6017 5969 int* sexternaldoflist=NULL; 6018 bool s ymmetric=true;5970 bool square=true; 6019 5971 6020 5972 /*retrieve some parameters: */ … … 6037 5989 } 6038 5990 6039 /*Use s ymmetricconstructor for ElementMatrix: */6040 if(!kff) Ke=new ElementMatrix(gsize,s ymmetric,gexternaldoflist);6041 else Ke=new ElementMatrix(gsize,s ymmetric,finternaldoflist,fexternaldoflist,fsize,sinternaldoflist,sexternaldoflist,ssize);5991 /*Use square constructor for ElementMatrix: */ 5992 if(!kff) Ke=new ElementMatrix(gsize,square,gexternaldoflist); 5993 else Ke=new ElementMatrix(gsize,square,finternaldoflist,fexternaldoflist,fsize,sinternaldoflist,sexternaldoflist,ssize); 6042 5994 6043 5995 /*Free ressources and return:*/ … … 6081 6033 } 6082 6034 6083 /* Use symmetricconstructor for ElementVector: */6035 /*constructor for ElementVector: */ 6084 6036 if(!kff)pe=new ElementVector(gsize,gexternaldoflist); 6085 6037 else pe=new ElementVector(gsize,finternaldoflist,fexternaldoflist,fsize); -
issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp
r5788 r5821 25 25 this->nrows=0; 26 26 this->ncols=0; 27 this->symmetric=false; 27 this->values=NULL; 28 this->square=false; 29 this->kff=false; 30 28 31 this->row_fsize=0; 29 this->values=NULL;30 32 this->row_finternaldoflist=NULL; 31 33 this->row_fexternaldoflist=NULL; … … 33 35 this->row_sinternaldoflist=NULL; 34 36 this->row_sexternaldoflist=NULL; 37 35 38 this->col_fsize=0; 36 39 this->col_finternaldoflist=NULL; … … 42 45 } 43 46 /*}}}*/ 44 /*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool s ymmetric,int* in_gexternaldoflist){{{1*/45 ElementMatrix::ElementMatrix(int gsize,bool s ymmetric,int* in_gexternaldoflist){46 47 if(s ymmetric=false)ISSMERROR(" calling symmetric constructor with false symmetricflag!");47 /*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool square,int* in_gexternaldoflist){{{1*/ 48 ElementMatrix::ElementMatrix(int gsize,bool square,int* in_gexternaldoflist){ 49 50 if(square=false)ISSMERROR(" calling square constructor with false square flag!"); 48 51 49 52 this->nrows=gsize; 50 53 this->ncols=gsize; 51 this->s ymmetric=true;54 this->square=true; 52 55 this->kff=false; 53 56 … … 59 62 memcpy(this->gexternaldoflist,in_gexternaldoflist,this->nrows*sizeof(int)); 60 63 } 61 else this->gexternaldoflist=NULL; 64 else{ 65 this->gexternaldoflist=NULL; 66 } 62 67 63 68 /*don't do rows and cols, because kff is false:*/ … … 78 83 } 79 84 /*}}}*/ 80 /*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool s ymmetric,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize){{{1*/81 ElementMatrix::ElementMatrix(int gsize,bool s ymmetric,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize){82 83 if(s ymmetric=false)ISSMERROR(" calling symmetric constructor with false symmetricflag!");85 /*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool square,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize){{{1*/ 86 ElementMatrix::ElementMatrix(int gsize,bool square,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize){ 87 88 if(square=false)ISSMERROR(" calling square constructor with false square flag!"); 84 89 85 90 this->nrows=gsize; 86 91 this->ncols=gsize; 87 this->s ymmetric=true;92 this->square=true; 88 93 this->kff=true; 89 94 … … 149 154 void ElementMatrix::AddValues(double* Ke_gg){ 150 155 151 int i,j;152 153 156 if(Ke_gg){ 154 for (i =0;i<this->nrows;i++){155 for( j=0;j<this->ncols;j++){157 for (int i=0;i<this->nrows;i++){ 158 for(int j=0;j<this->ncols;j++){ 156 159 *(this->values+this->ncols*i+j)+=*(Ke_gg+this->ncols*i+j); 157 160 } … … 166 169 double* internalvalues=NULL; 167 170 168 if(this->s ymmetric){171 if(this->square){ 169 172 /*only use row dofs to add values into global matrices: */ 170 173 … … 207 210 } 208 211 else{ 209 ISSMERROR(" unsymmetricAddToGlobal routine not support yet!");212 ISSMERROR(" non square matrix AddToGlobal routine not support yet!"); 210 213 } 211 214 … … 219 222 printf(" nrows: %i\n",nrows); 220 223 printf(" ncols: %i\n",ncols); 221 printf(" s ymmetric: %s\n",symmetric?"true":"false");224 printf(" square: %s\n",square?"true":"false"); 222 225 printf(" kff: %s\n",kff?"true":"false"); 223 226 … … 225 228 for(i=0;i<nrows;i++){ 226 229 printf(" %i: ",i); 227 for(j=0;j<ncols;j++){ 228 printf("%g ",*(values+ncols*i+j)); 229 } 230 for(j=0;j<ncols;j++) printf("%g ",*(values+ncols*i+j)); 230 231 printf("\n"); 231 232 } … … 246 247 if(row_sexternaldoflist)for(i=0;i<row_ssize;i++)printf("%i ",row_sexternaldoflist[i]); printf("\n"); 247 248 248 if(!symmetric){ 249 249 if(!square){ 250 250 printf(" col_fsize: %i\n",col_fsize); 251 251 printf(" col_finternaldoflist (%p): ",col_finternaldoflist); -
issm/trunk/src/c/objects/Numerics/ElementMatrix.h
r5772 r5821 22 22 int ncols; 23 23 double* values; 24 bool s ymmetric;24 bool square; 25 25 bool kff; 26 26 … … 50 50 /*ElementMatrix constructors, destructors {{{1*/ 51 51 ElementMatrix(); 52 ElementMatrix(int gsize,bool s ymmetric,int* gexternaldoflist);53 ElementMatrix(int gsize,bool s ymmetric,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize);52 ElementMatrix(int gsize,bool square,int* gexternaldoflist); 53 ElementMatrix(int gsize,bool square,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize); 54 54 ~ElementMatrix(); 55 55 /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.