Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 2704) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 2705) @@ -24,12 +24,15 @@ //#define _DEBUGGAUSS_ +/*FUNCTION Tria constructor {{{1*/ Tria::Tria(){ return; } - +/*}}}*/ +/*FUNCTION Tria destructor {{{1*/ Tria::~Tria(){ return; } - +/*}}}*/ +/*FUNCTION Tria creation {{{1*/ Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_numparid,int tria_node_ids[3],double tria_h[3],double tria_s[3],double tria_b[3],double tria_k[3],double tria_melting[3], double tria_accumulation[3],double tria_geothermalflux[3],int tria_friction_type,double tria_p,double tria_q,int tria_shelf, bool tria_onwater){ @@ -68,5 +71,6 @@ return; } - +/*}}}*/ +/*FUNCTION Echo {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Tria::Echo" @@ -104,4 +108,6 @@ return; } +/*}}}*/ +/*FUNCTION DeepEcho{{{1*/ #undef __FUNCT__ #define __FUNCT__ "Tria::DeepEcho" @@ -139,4 +145,6 @@ return; } +/*}}}*/ +/*FUNCTION Marshall {{{1*/ void Tria::Marshall(char** pmarshalled_dataset){ @@ -184,5 +192,6 @@ return; } - +/*}}}*/ +/*FUNCTION MarshallSize {{{1*/ int Tria::MarshallSize(){ return sizeof(id) @@ -214,9 +223,11 @@ +sizeof(int); //sizeof(int) for enum type } - +/*}}}*/ +/*FUNCTION GetName {{{1*/ char* Tria::GetName(void){ return "tria"; } - +/*}}}*/ +/*FUNCTION Demarshall {{{1*/ void Tria::Demarshall(char** pmarshalled_dataset){ @@ -267,4 +278,6 @@ return; } +/*}}}*/ +/*FUNCTION Enum {{{1*/ int Tria::Enum(void){ @@ -272,12 +285,15 @@ } +/*}}}*/ +/*FUNCTION GetId {{{1*/ int Tria::GetId(){ return id; } - +/*}}}*/ +/*FUNCTION MyRank {{{1*/ int Tria::MyRank(void){ extern int my_rank; return my_rank; } - - +/*}}}*/ +/*FUNCTION Configure {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Tria::Configure" @@ -308,5 +324,6 @@ } - +/*}}}*/ +/*FUNCTION CreateKMatrix {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Tria::CreateKMatrix" @@ -343,6 +360,6 @@ } - - +/*}}}*/ +/*FUNCTION CreateKMatrixDiagnosticHoriz {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Tria::CreateKMatrixDiagnosticHoriz" @@ -562,4 +579,6 @@ } +/*}}}*/ +/*FUNCTION CreateKMatrixPrognostic {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Tria::CreateKMatrixPrognostic" @@ -780,8 +799,9 @@ } - -#undef __FUNCT__ -#define __FUNCT__ "Tria::CreatePVectorPrognostic" -void Tria::CreatePVectorPrognostic(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){ +/*}}}*/ +/*FUNCTION CreateKMatrixDiagnosticHorizFriction {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Tria::CreateKMatrixDiagnosticHorizFriction" +void Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ @@ -791,6 +811,5 @@ /* node data: */ const int numgrids=3; - const int NDOF1=1; - const int numdof=NDOF1*numgrids; + const int numdof=2*numgrids; double xyz_list[numgrids][3]; int doflist[numdof]; @@ -806,9 +825,466 @@ double gauss_l1l2l3[3]; + /* matrices: */ + double L[2][numdof]; + double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag + double DL_scalar; + + /* local element matrices: */ + double Ke_gg[numdof][numdof]; //local element stiffness matrix + double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag + + double Jdet; + + /*slope: */ + double slope[2]={0.0,0.0}; + double slope_magnitude; + + /*input parameters for structural analysis (diagnostic): */ + double vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; + int dofs[2]={0,1}; + + /*friction: */ + double alpha2_list[numgrids]={0.0,0.0,0.0}; + double alpha2; + + double MAXSLOPE=.06; // 6 % + double MOUNTAINKEXPONENT=10; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Set Ke_gg to 0: */ + for(i=0;iRecover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + + /*Build alpha2_list used by drag stiffness matrix*/ + Friction* friction=NewFriction(); + + /*Initialize all fields: */ + friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char)); + strcpy(friction->element_type,"2d"); + + friction->gravity=matpar->GetG(); + friction->rho_ice=matpar->GetRhoIce(); + friction->rho_water=matpar->GetRhoWater(); + friction->K=&k[0]; + friction->bed=&b[0]; + friction->thickness=&h[0]; + friction->velocities=&vxvy_list[0][0]; + friction->p=p; + friction->q=q; + + /*Compute alpha2_list: */ + FrictionGetAlpha2(&alpha2_list[0],friction); + + /*Erase friction object: */ + DeleteFriction(&friction); + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && id==ELID){ + printf(" alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]); + } + #endif + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && id==ELID){ + printf(" gaussian points: \n"); + for(i=0;i 6% for this element, it means we are on a mountain. In this particular case, + //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ + GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3); + slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); + + if (slope_magnitude>MAXSLOPE){ + alpha2_list[0]=pow((double)10,MOUNTAINKEXPONENT); + alpha2_list[1]=pow((double)10,MOUNTAINKEXPONENT); + alpha2_list[2]=pow((double)10,MOUNTAINKEXPONENT); + } + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /*Get L matrix: */ + GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); + + /*Now, take care of the basal friction if there is any: */ + GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3); + + DL_scalar=alpha2*gauss_weight*Jdet; + for (i=0;i<2;i++){ + DL[i][i]=DL_scalar; + } + + /* Do the triple producte tL*D*L: */ + TripleMultiply( &L[0][0],2,numdof,1, + &DL[0][0],2,2,0, + &L[0][0],2,numdof,0, + &Ke_gg_gaussian[0][0],0); + + for( i=0; iGetG()); + printf("rho_ice %g\n",matpar->GetRhoIce()); + printf("thickness [%g,%g,%g]\n",h[0],h[1],h[2]); + printf("surface[%g,%g,%g]\n",s[0],s[1],s[2]); + printf("bed[%g,%g,%g]\n",b[0],b[1],b[2]); + printf("drag [%g,%g,%g]\n",k[0],k[1],k[2]); + } + #endif + + + /* Get gaussian points and weights: */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/ + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && id==ELID){ + printf(" gaussian points: \n"); + for(i=0;iGetRhoIce()*matpar->GetG()*thickness; + + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && id==ELID){ + printf(" gaussian %i\n",ig); + printf(" thickness %g\n",thickness); + printf(" slope(%g,%g)\n",slope[0],slope[1]); + printf(" Jdet %g\n",Jdet); + printf(" gaussweigth %g\n",gauss_weight); + printf(" l1l2l3 (%g,%g,%g)\n",l1l2l3[0],l1l2l3[1],l1l2l3[2]); + if(friction_type==1)printf(" plastic_stress(%g)\n",plastic_stress); + } + #endif + + /*Build pe_g_gaussian vector: */ + if(friction_type==1){ + for (i=0;iterms\n",ig); + for( i=0; inrows; i++){ + printf("%g ",*(pe_g->terms+i)); + } + printf("\n"); + printf(" pe_g->row_indices\n",ig); + for( i=0; inrows; i++){ + printf("%i ",*(pe_g->row_indices+i)); + } + } + #endif + + /*Add pe_g to global vector pg: */ + VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + +} +/*}}}*/ +/*FUNCTION CreatePVectorPrognostic {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Tria::CreatePVectorPrognostic" +void Tria::CreatePVectorPrognostic(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){ + + + /* local declarations */ + int i,j; + + /* node data: */ + const int numgrids=3; + const int NDOF1=1; + const int numdof=NDOF1*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdof]; + int numberofdofspernode; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + double gauss_weight; + double gauss_l1l2l3[3]; + /* matrix */ double pe_g[numgrids]={0.0}; double L[numgrids]; double Jdettria; - + /*input parameters for structural analysis (diagnostic): */ double accumulation_list[numgrids]={0.0}; @@ -862,8 +1338,8 @@ GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3); GetParameterValue(&thickness_g, &thickness_list[0],gauss_l1l2l3); - + /* Add value into pe_g: */ for( i=0; iRecover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); - - /*Build alpha2_list used by drag stiffness matrix*/ - Friction* friction=NewFriction(); - - /*Initialize all fields: */ - friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char)); - strcpy(friction->element_type,"2d"); - - friction->gravity=matpar->GetG(); - friction->rho_ice=matpar->GetRhoIce(); - friction->rho_water=matpar->GetRhoWater(); - friction->K=&k[0]; - friction->bed=&b[0]; - friction->thickness=&h[0]; - friction->velocities=&vxvy_list[0][0]; - friction->p=p; - friction->q=q; - - /*Compute alpha2_list: */ - FrictionGetAlpha2(&alpha2_list[0],friction); - - /*Erase friction object: */ - DeleteFriction(&friction); - - #ifdef _DEBUGELEMENTS_ - if(my_rank==RANK && id==ELID){ - printf(" alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]); - } - #endif - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - #ifdef _DEBUGELEMENTS_ - if(my_rank==RANK && id==ELID){ - printf(" gaussian points: \n"); - for(i=0;i 6% for this element, it means we are on a mountain. In this particular case, - //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ - GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3); - slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); - - if (slope_magnitude>MAXSLOPE){ - alpha2_list[0]=pow((double)10,MOUNTAINKEXPONENT); - alpha2_list[1]=pow((double)10,MOUNTAINKEXPONENT); - alpha2_list[2]=pow((double)10,MOUNTAINKEXPONENT); - } - + GetParameterDerivativeValue(&slope[0], ¶m[0],&xyz_list[0][0], gauss_l1l2l3); + /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /*Get L matrix: */ - GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); - - /*Now, take care of the basal friction if there is any: */ - GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3); - - DL_scalar=alpha2*gauss_weight*Jdet; - for (i=0;i<2;i++){ - DL[i][i]=DL_scalar; - } - - /* Do the triple producte tL*D*L: */ - TripleMultiply( &L[0][0],2,numdof,1, - &DL[0][0],2,2,0, - &L[0][0],2,numdof,0, - &Ke_gg_gaussian[0][0],0); - - for( i=0; iGetG()); - printf("rho_ice %g\n",matpar->GetRhoIce()); - printf("thickness [%g,%g,%g]\n",h[0],h[1],h[2]); - printf("surface[%g,%g,%g]\n",s[0],s[1],s[2]); - printf("bed[%g,%g,%g]\n",b[0],b[1],b[2]); - printf("drag [%g,%g,%g]\n",k[0],k[1],k[2]); - } - #endif - - - /* Get gaussian points and weights: */ - GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/ - - #ifdef _DEBUGELEMENTS_ - if(my_rank==RANK && id==ELID){ - printf(" gaussian points: \n"); - for(i=0;iGetRhoIce()*matpar->GetG()*thickness; - - - #ifdef _DEBUGELEMENTS_ - if(my_rank==RANK && id==ELID){ - printf(" gaussian %i\n",ig); - printf(" thickness %g\n",thickness); - printf(" slope(%g,%g)\n",slope[0],slope[1]); - printf(" Jdet %g\n",Jdet); - printf(" gaussweigth %g\n",gauss_weight); - printf(" l1l2l3 (%g,%g,%g)\n",l1l2l3[0],l1l2l3[1],l1l2l3[2]); - if(friction_type==1)printf(" plastic_stress(%g)\n",plastic_stress); - } - #endif - - /*Build pe_g_gaussian vector: */ - if(friction_type==1){ - for (i=0;iterms\n",ig); - for( i=0; inrows; i++){ - printf("%g ",*(pe_g->terms+i)); - } - printf("\n"); - printf(" pe_g->row_indices\n",ig); - for( i=0; inrows; i++){ - printf("%i ",*(pe_g->row_indices+i)); - } - } - #endif - - /*Add pe_g to global vector pg: */ - VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); - - cleanup_and_return: - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); - -} - -#undef __FUNCT__ -#define __FUNCT__ "Tria::CreatePVectorSlopeCompute" - -void Tria::CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){ - - int i,j; - - /* node data: */ - const int numgrids=3; - const int NDOF1=1; - const int numdof=NDOF1*numgrids; - double xyz_list[numgrids][3]; - int doflist[numdof]; - int numberofdofspernode; - - /* gaussian points: */ - int num_gauss,ig; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* gauss_weights = NULL; - double gauss_weight; - double gauss_l1l2l3[3]; - - /* Jacobian: */ - double Jdet; - - /*nodal functions: */ - double l1l2l3[3]; - - /*element vector at the gaussian points: */ - double pe_g[numdof]; - double pe_g_gaussian[numdof]; - double param[3]; - double slope[2]; - - /* Get node coordinates and dof list: */ - GetElementNodeData( &xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* Set pe_g to 0: */ - for(i=0;i