Changeset 5654
- Timestamp:
- 09/02/10 09:47:07 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Penta.cpp
r5651 r5654 2607 2607 void Penta::CreateKMatrixDiagnosticMacAyeal( Mat Kgg){ 2608 2608 2609 /* local declarations */ 2610 int i,j; 2611 2612 /* node data: */ 2613 const int numgrids3d=6; 2614 const int numgrids2d=3; 2615 const int numdof2d=2*numgrids2d; 2616 double xyz_list[numgrids3d][3]; 2617 int* doflist=NULL; 2618 2619 /* 3d gaussian points: */ 2620 int ig; 2621 GaussPenta *gauss=NULL; 2622 GaussTria *gauss_tria=NULL; 2623 2624 /* material data: */ 2625 double viscosity; //viscosity 2626 double oldviscosity; //viscosity 2627 double newviscosity; //viscosity 2628 2629 /* strain rate: */ 2630 double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ 2631 double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ 2632 2633 /* matrices: */ 2634 double B[3][numdof2d]; 2635 double Bprime[3][numdof2d]; 2636 double L[2][numdof2d]; 2637 double D[3][3]={0.0}; // material matrix, simple scalar matrix. 2638 double D_scalar; 2639 double DL[2][2]={0.0}; //for basal drag 2640 double DL_scalar; 2641 2642 /* local element matrices: */ 2643 double Ke_gg[numdof2d][numdof2d]={0.0}; //local element stiffness matrix 2644 double Ke_gg_gaussian[numdof2d][numdof2d]; //stiffness matrix evaluated at the gaussian point. 2645 double Jdet; 2646 2647 /*slope: */ 2648 double slope[2]={0.0}; 2649 double slope_magnitude; 2650 2651 /*friction: */ 2652 double alpha2_list[3]; 2653 double alpha2; 2654 2655 double MAXSLOPE=.06; // 6 % 2656 double MOUNTAINKEXPONENT=10; 2657 2658 /*parameters: */ 2659 double viscosity_overshoot; 2660 2609 2661 /*Collapsed formulation: */ 2610 2662 Tria* tria=NULL; 2663 Penta* pentabase=NULL; 2611 2664 2612 2665 /*inputs: */ 2613 2666 bool onwater; 2614 2667 bool onbed; 2615 2668 bool shelf; 2669 int approximation; 2670 Input* vx_input=NULL; 2671 Input* vy_input=NULL; 2672 Input* vxold_input=NULL; 2673 Input* vyold_input=NULL; 2674 2675 inputs->GetParameterValue(&approximation,ApproximationEnum); 2616 2676 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 2617 2677 inputs->GetParameterValue(&onbed,ElementOnBedEnum); 2678 inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); 2618 2679 2619 2680 /*If on water, skip stiffness: */ 2620 if(onwater)return; 2621 2622 /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the 2623 bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build 2624 the stiffness matrix. */ 2625 2626 if (onbed==0){ 2627 /*This element should be collapsed, but this element is not on the bedrock, therefore all its 2628 * dofs have already been frozen! Do nothing: */ 2629 return; 2630 } 2631 else if (onbed==1){ 2632 2633 /*This element should be collapsed into a tria element at its base. Create this tria element, 2634 *and use its CreateKMatrix functionality to fill the global stiffness matrix: */ 2635 2636 /*Depth Averaging B*/ 2637 this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum); 2638 2639 /*Call Tria function*/ 2640 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 2641 tria->CreateKMatrix(Kgg); 2642 delete tria->matice; delete tria; 2643 2644 /*Delete B averaged*/ 2645 this->matice->inputs->DeleteInput(RheologyBbarEnum); 2646 2647 return; 2681 if(approximation==MacAyealApproximationEnum){ 2682 if(onwater)return; 2683 2684 /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the 2685 bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build 2686 the stiffness matrix. */ 2687 2688 if (onbed==0){ 2689 /*This element should be collapsed, but this element is not on the bedrock, therefore all its 2690 * dofs have already been frozen! Do nothing: */ 2691 return; 2692 } 2693 else if (onbed==1){ 2694 2695 /*This element should be collapsed into a tria element at its base. Create this tria element, 2696 *and use its CreateKMatrix functionality to fill the global stiffness matrix: */ 2697 2698 /*Depth Averaging B*/ 2699 this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum); 2700 2701 /*Call Tria function*/ 2702 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 2703 tria->CreateKMatrix(Kgg); 2704 delete tria->matice; delete tria; 2705 2706 /*Delete B averaged*/ 2707 this->matice->inputs->DeleteInput(RheologyBbarEnum); 2708 2709 return; 2710 } 2711 } 2712 else if(approximation==MacAyealPattynApproximationEnum){ 2713 /*retrieve some parameters: */ 2714 this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); 2715 2716 /*If on water, skip stiffness: */ 2717 if(onwater)return; 2718 2719 /*Find penta on bed as this is a macayeal elements: */ 2720 pentabase=GetBasalElement(); 2721 tria=pentabase->SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. 2722 2723 /* Get node coordinates and dof list: */ 2724 GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids3d); 2725 tria->GetDofList(&doflist,MacAyealApproximationEnum); //Pattyn dof list 2726 2727 /*Retrieve all inputs we will be needing: */ 2728 vx_input=inputs->GetInput(VxEnum); 2729 vy_input=inputs->GetInput(VyEnum); 2730 vxold_input=inputs->GetInput(VxOldEnum); 2731 vyold_input=inputs->GetInput(VyOldEnum); 2732 2733 /* Start looping on the number of gaussian points: */ 2734 gauss=new GaussPenta(5,5); 2735 gauss_tria=new GaussTria(); 2736 for (ig=gauss->begin();ig<gauss->end();ig++){ 2737 2738 gauss->GaussPoint(ig); 2739 gauss->SynchronizeGaussTria(gauss_tria); 2740 2741 /*Get strain rate from velocity: */ 2742 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); 2743 this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); 2744 2745 /*Get viscosity: */ 2746 matice->GetViscosity3d(&viscosity, &epsilon[0]); 2747 matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); 2748 2749 /*Get B and Bprime matrices: */ 2750 tria->GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss_tria); 2751 tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria); 2752 2753 /* Get Jacobian determinant: */ 2754 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); 2755 2756 /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant 2757 onto this scalar matrix, so that we win some computational time: */ 2758 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); 2759 D_scalar=2*newviscosity*gauss->weight*Jdet; 2760 for (i=0;i<3;i++) D[i][i]=D_scalar; 2761 2762 /* Do the triple product tB*D*Bprime: */ 2763 TripleMultiply( &B[0][0],3,numdof2d,1, 2764 &D[0][0],3,3,0, 2765 &Bprime[0][0],3,numdof2d,0, 2766 &Ke_gg_gaussian[0][0],0); 2767 2768 /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ 2769 for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j]; 2770 } 2771 2772 /*Add Ke_gg to global matrix Kgg: */ 2773 MatSetValues(Kgg,numdof2d,doflist,numdof2d,doflist,(const double*)Ke_gg,ADD_VALUES); 2774 2775 //Deal with 2d friction at the bedrock interface 2776 if((onbed && !shelf)){ 2777 /*Build a tria element using the 3 grids of the base of the penta. Then use 2778 * the tria functionality to build a friction stiffness matrix on these 3 2779 * grids: */ 2780 2781 tria->CreateKMatrixDiagnosticMacAyealFriction(Kgg); 2782 } 2783 2784 /*Clean up and return*/ 2785 delete tria->matice; 2786 delete tria; 2787 delete gauss_tria; 2788 delete gauss; 2789 xfree((void**)&doflist); 2648 2790 } 2649 2791 }
Note:
See TracChangeset
for help on using the changeset viewer.