Changeset 5840
- Timestamp:
- 09/16/10 10:14:40 (15 years ago)
- Location:
- issm/trunk/src/c/objects/Elements
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Tria.cpp
r5839 r5840 695 695 void Tria::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){ 696 696 697 /*retrive parameters: */ 697 /*retreive parameters: */ 698 ElementMatrix* Ke=NULL; 698 699 int analysis_type; 699 700 parameters->FindParam(&analysis_type,AnalysisTypeEnum); … … 704 705 705 706 /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ 706 if (analysis_type==DiagnosticHorizAnalysisEnum || analysis_type==AdjointHorizAnalysisEnum){ 707 ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal(); 707 switch(analysis_type){ 708 case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum: 709 Ke=CreateKMatrixDiagnosticMacAyeal(); 710 break; 711 case DiagnosticHutterAnalysisEnum: 712 Ke=CreateKMatrixDiagnosticHutter(); 713 break; 714 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: 715 Ke=CreateKMatrixSlope(); 716 break; 717 case PrognosticAnalysisEnum: 718 Ke=CreateKMatrixPrognostic(); 719 break; 720 case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum: 721 Ke=CreateKMatrixBalancedthickness(); 722 break; 723 case BalancedvelocitiesAnalysisEnum: 724 Ke=CreateKMatrixBalancedvelocities(); 725 break; 726 default: 727 ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); 728 } 729 730 /*Add to global matrix*/ 731 if(Ke){ 708 732 Ke->AddToGlobal(Kgg,Kff,Kfs); 709 733 delete Ke; 710 734 } 711 else if (analysis_type==DiagnosticHutterAnalysisEnum){712 CreateKMatrixDiagnosticHutter( Kgg);713 }714 else if (analysis_type==BedSlopeXAnalysisEnum || analysis_type==SurfaceSlopeXAnalysisEnum || analysis_type==BedSlopeYAnalysisEnum || analysis_type==SurfaceSlopeYAnalysisEnum){715 CreateKMatrixSlope( Kgg);716 }717 else if (analysis_type==PrognosticAnalysisEnum){718 ElementMatrix* Ke=CreateKMatrixPrognostic();719 Ke->AddToGlobal(Kgg,Kff,Kfs);720 delete Ke;721 }722 else if (analysis_type==BalancedthicknessAnalysisEnum || analysis_type==AdjointBalancedthicknessAnalysisEnum){723 ElementMatrix* Ke=CreateKMatrixBalancedthickness();724 Ke->AddToGlobal(Kgg,Kff,Kfs);725 delete Ke;726 }727 else if (analysis_type==BalancedvelocitiesAnalysisEnum){728 CreateKMatrixBalancedvelocities( Kgg);729 }730 else{731 ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));732 }733 734 735 } 735 736 /*}}}*/ … … 2494 2495 &Ke_gg_thickness2[0][0],0); 2495 2496 2496 for( i=0; i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_thickness1[i][j];2497 for( i=0; i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_thickness2[i][j];2497 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_thickness1[i][j]; 2498 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_thickness2[i][j]; 2498 2499 2499 2500 if(artdiff){ … … 2506 2507 &Ke_gg_gaussian[0][0],0); 2507 2508 2508 for( i=0; i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j];2509 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j]; 2509 2510 } 2510 2511 } … … 2578 2579 /*}}}*/ 2579 2580 /*FUNCTION Tria::CreateKMatrixBalancedvelocities {{{1*/ 2580 void Tria::CreateKMatrixBalancedvelocities(Mat Kgg){2581 ElementMatrix* Tria::CreateKMatrixBalancedvelocities(void){ 2581 2582 2582 2583 /* local declarations */ 2584 const int numdof=NDOF1*NUMVERTICES; 2583 2585 int i,j; 2584 2585 /* node data: */2586 const int numdof=NDOF1*NUMVERTICES;2587 2586 double xyz_list[NUMVERTICES][3]; 2588 int* doflist=NULL;2589 2590 /* gaussian points: */2591 2587 int ig; 2592 2588 GaussTria *gauss=NULL; 2593 2594 /* matrices: */2595 2589 double L[NUMVERTICES]; 2596 2590 double B[2][NUMVERTICES]; … … 2604 2598 double Ke_gg_velocities2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. 2605 2599 double Jdettria; 2606 2607 /*input parameters for structural analysis (diagnostic): */2608 2600 double surface_normal[3]; 2609 2601 double surface_list[3]; … … 2617 2609 double KDL[2][2]={0.0}; 2618 2610 int dim; 2619 2620 /*inputs: */ 2611 bool artdiff; 2612 2613 /*Initialize Element matrix and return if necessary*/ 2614 if(IsOnWater()) return NULL; 2615 ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum); 2616 2617 /*Retrieve all inputs and parameters*/ 2618 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 2619 this->parameters->FindParam(&artdiff,ArtDiffEnum); 2620 this->parameters->FindParam(&dim,DimEnum); 2621 Input* surface_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input); 2621 2622 Input* vxaverage_input=NULL; 2622 2623 Input* vyaverage_input=NULL; 2623 bool artdiff;2624 2625 /*retrieve some parameters: */2626 this->parameters->FindParam(&artdiff,ArtDiffEnum);2627 this->parameters->FindParam(&dim,DimEnum);2628 2629 /* Get node coordinates and dof list: */2630 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);2631 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);2632 2633 /*Retrieve all inputs we will be needed*/2634 Input* surface_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input);2635 2624 if(dim==2){ 2636 2625 vxaverage_input=inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input); … … 2664 2653 //Create Artificial diffusivity once for all if requested 2665 2654 if(artdiff){ 2666 //Get the Jacobian determinant2667 2655 gauss=new GaussTria(); 2668 2656 gauss->GaussCenter(); … … 2670 2658 delete gauss; 2671 2659 2672 //Build K matrix (artificial diffusivity matrix)2673 2660 vxaverage_input->GetParameterAverage(&v_gauss[0]); 2674 2661 vyaverage_input->GetParameterAverage(&v_gauss[1]); … … 2684 2671 gauss->GaussPoint(ig); 2685 2672 2686 /* Get Jacobian determinant: */2687 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);2688 2689 /*Get B and B prime matrix: */2690 2673 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss); 2691 2674 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss); 2692 2693 //Get vx, vy and their derivatives at gauss point 2675 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); 2676 2694 2677 vxaverage_input->GetParameterValue(&vx,gauss); 2695 2678 vyaverage_input->GetParameterValue(&vy,gauss); 2696 2697 2679 vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); 2698 2680 vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); 2699 2700 2681 dvxdx=dvx[0]; 2701 2682 dvydy=dvy[1]; … … 2706 2687 DLprime[1][1]=DL_scalar*ny; 2707 2688 2708 //Do the triple product tL*D*L.2709 //Ke_gg_velocities=B'*DLprime*Bprime;2710 2689 TripleMultiply( &B[0][0],2,numdof,1, 2711 2690 &DLprime[0][0],2,2,0, … … 2713 2692 &Ke_gg_velocities2[0][0],0); 2714 2693 2715 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ 2716 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_velocities2[i][j]; 2694 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_velocities2[i][j]; 2717 2695 2718 2696 if(artdiff){ 2719 2720 /* Compute artificial diffusivity */2721 2697 KDL[0][0]=DL_scalar*K[0][0]; 2722 2698 KDL[1][1]=DL_scalar*K[1][1]; … … 2727 2703 &Ke_gg_gaussian[0][0],0); 2728 2704 2729 /* Add artificial diffusivity matrix */ 2730 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j]; 2705 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j]; 2731 2706 2732 2707 } 2733 2708 2734 2709 } 2735 2736 /*Add Ke_gg to global matrix Kgg: */2737 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);2738 2710 2739 2711 /*Clean up and return*/ 2740 2712 delete gauss; 2741 xfree((void**)&doflist);2713 return Ke; 2742 2714 } 2743 2715 /*}}}*/ … … 3139 3111 /*}}}*/ 3140 3112 /*FUNCTION Tria::CreateKMatrixDiagnosticHutter{{{1*/ 3141 void Tria::CreateKMatrixDiagnosticHutter(Mat Kgg){ 3142 3143 /*Collapsed formulation: */ 3144 int i; 3145 int connectivity; 3146 const int numdofs=NUMVERTICES*NDOF2; 3147 int* doflist=NULL; 3148 double Ke_gg[numdofs][numdofs]={0.0}; 3149 3150 /*If on water, skip: */ 3151 if(IsOnWater())return; 3152 3153 GetDofList(&doflist,NoneApproximationEnum,GsetEnum); 3154 3155 /*Spawn 3 sing elements: */ 3156 for(i=0;i<3;i++){ 3113 ElementMatrix* Tria::CreateKMatrixDiagnosticHutter(void){ 3114 3115 /*Intermediaries*/ 3116 const int numdof=NUMVERTICES*NDOF2; 3117 int i,connectivity; 3118 3119 /*Initialize Element matrix and return if necessary*/ 3120 if(IsOnWater()) return NULL; 3121 ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum); 3122 3123 /*Create Element matrix*/ 3124 for(i=0;i<NUMVERTICES;i++){ 3157 3125 connectivity=nodes[i]->GetConnectivity(); 3158 Ke_gg[2*i][2*i]=1/(double)connectivity; 3159 Ke_gg[2*i+1][2*i+1]=1/(double)connectivity; 3160 } 3161 3162 MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES); 3163 3164 /*Free ressources:*/ 3165 xfree((void**)&doflist); 3126 Ke->values[(2*i)*numdof +(2*i) ]=1/(double)connectivity; 3127 Ke->values[(2*i+1)*numdof+(2*i+1)]=1/(double)connectivity; 3128 } 3129 3130 /*Clean up and return*/ 3131 return Ke; 3166 3132 } 3167 3133 /*}}}*/ … … 3570 3536 /*}}}*/ 3571 3537 /*FUNCTION Tria::CreateKMatrixSlope {{{1*/ 3572 void Tria::CreateKMatrixSlope(Mat Kgg){3538 ElementMatrix* Tria::CreateKMatrixSlope(void){ 3573 3539 3574 3540 /*constants: */ 3575 const int numnodes=3; 3576 const int numdof=NDOF1*numnodes; 3541 const int numdof=NDOF1*NUMVERTICES; 3577 3542 3578 3543 /* Intermediaries */ 3579 3544 int i,j,ig; 3580 3545 double DL_scalar,Jdet; 3581 double xyz_list[ numnodes][3];3546 double xyz_list[NUMVERTICES][3]; 3582 3547 double L[1][3]; 3583 double Ke[numdof][numdof] = {0.0};3584 3548 double Ke_g[numdof][numdof]; 3585 3549 GaussTria *gauss = NULL; 3586 int *doflist = NULL; 3587 3588 GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes); 3589 GetDofList(&doflist,NoneApproximationEnum,GsetEnum); 3550 3551 /*Initialize Element matrix and return if necessary*/ 3552 if(IsOnWater()) return NULL; 3553 ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum); 3554 3555 GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES); 3590 3556 3591 3557 /* Start looping on the number of gaussian points: */ … … 3605 3571 &Ke_g[0][0],0); 3606 3572 3607 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke[i][j]+=Ke_g[i][j]; 3608 } 3609 3610 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES); 3611 3612 /*Clean up*/ 3573 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j]; 3574 } 3575 3576 /*Clean up and return*/ 3613 3577 delete gauss; 3614 xfree((void**)&doflist);3578 return Ke; 3615 3579 } 3616 3580 /*}}}*/ -
issm/trunk/src/c/objects/Elements/Tria.h
r5839 r5840 124 124 ElementMatrix* CreateKMatrixBalancedthickness_DG(void); 125 125 ElementMatrix* CreateKMatrixBalancedthickness_CG(void); 126 void CreateKMatrixBalancedvelocities(Mat Kgg);126 ElementMatrix* CreateKMatrixBalancedvelocities(void); 127 127 ElementMatrix* CreateKMatrixDiagnosticMacAyeal(void); 128 128 ElementMatrix* CreateKMatrixDiagnosticMacAyealViscous(void); … … 131 131 void CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg); 132 132 void CreateKMatrixDiagnosticPattynFriction(Mat Kgg); 133 void CreateKMatrixDiagnosticHutter(Mat Kgg);133 ElementMatrix* CreateKMatrixDiagnosticHutter(void); 134 134 void CreateKMatrixDiagnosticSurfaceVert(Mat Kgg); 135 135 void CreateKMatrixMelting(Mat Kgg); … … 137 137 ElementMatrix* CreateKMatrixPrognostic_CG(void); 138 138 ElementMatrix* CreateKMatrixPrognostic_DG(void); 139 void CreateKMatrixSlope(Mat Kgg);139 ElementMatrix* CreateKMatrixSlope(void); 140 140 void CreateKMatrixThermal(Mat Kgg); 141 141 void CreatePVectorBalancedthickness_DG(Vec pg);
Note:
See TracChangeset
for help on using the changeset viewer.