![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 5821 for issm/trunk/src/c/objects/Elements/Tria.cpp
- Timestamp:
- 09/15/10 10:18:38 (15 years ago)
- File:
-
- 1 edited
-
issm/trunk/src/c/objects/Elements/Tria.cpp (modified) (23 diffs)
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);
Note:
See TracChangeset
for help on using the changeset viewer.
