Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp (revision 393) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp (revision 394) @@ -91,5 +91,5 @@ /*Free data: */ xfree((void**)&gridonstokes); - + /*Assign output pointer: */ *pconstraints=constraints; Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (revision 393) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (revision 394) @@ -514,5 +514,5 @@ xfree((void**)&model->gridonstokes); xfree((void**)&model->borderstokes); - + cleanup_and_return: Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp (revision 393) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp (revision 394) @@ -51,4 +51,5 @@ int numberofsegs_diag_stokes; int count; + /*Create loads: */ @@ -134,15 +135,16 @@ #endif - if ((model->gridonbed[i]) && (model->gridonicesheet[i]) && (model->gridonstokes[i])) + if ((model->gridonbed[i]) && (model->gridonicesheet[i]) && (model->gridonstokes[i])){ - pengrid_id=count+1; //matlab indexing - pengrid_node_id=i+1; - pengrid_dof=1; - pengrid_penalty_offset=model->penalty_offset; + pengrid_id=count+1; //matlab indexing + pengrid_node_id=i+1; + pengrid_dof=1; + pengrid_penalty_offset=model->penalty_offset; - pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate); - - loads->AddObject(icefront); - count++; + pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate); + + loads->AddObject(icefront); + count++; + } #ifdef _PARALLEL_ } //if((model->my_grids[i]==1)) Index: /issm/trunk/src/c/ModelProcessorx/Model.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 393) +++ /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 394) @@ -307,4 +307,5 @@ ModelFetchData((void**)&model->ishutter,NULL,NULL,model_handle,"ishutter","Integer",NULL); ModelFetchData((void**)&model->ismacayealpattyn,NULL,NULL,model_handle,"ismacayealpattyn","Integer",NULL); + ModelFetchData((void**)&model->isstokes,NULL,NULL,model_handle,"isstokes","Integer",NULL); /*!Get drag_type, drag and p,q: */ Index: /issm/trunk/src/c/objects/Friction.cpp =================================================================== --- /issm/trunk/src/c/objects/Friction.cpp (revision 393) +++ /issm/trunk/src/c/objects/Friction.cpp (revision 394) @@ -35,4 +35,5 @@ friction->p=UNDEF; friction->q=UNDEF; + friction->analysis_type=UNDEF; return 1; Index: /issm/trunk/src/c/objects/Friction.h =================================================================== --- /issm/trunk/src/c/objects/Friction.h (revision 393) +++ /issm/trunk/src/c/objects/Friction.h (revision 394) @@ -19,4 +19,5 @@ double p; double q; + int analysis_type; }; Index: /issm/trunk/src/c/objects/Matice.cpp =================================================================== --- /issm/trunk/src/c/objects/Matice.cpp (revision 393) +++ /issm/trunk/src/c/objects/Matice.cpp (revision 394) @@ -306,4 +306,72 @@ } +#undef __FUNCT__ +#define __FUNCT__ "MatIce::GetViscosity3dStokes" +void Matice::GetViscosity3dStokes(double* pviscosity3d, double* epsilon){ + + /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: + * + * 2*B + * viscosity3d= ------------------------------------------------------------------- + * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] + * + * where mu is the viscotiy, B the flow law parameter , (u,v) the velocity + * vector, and n the flow law exponent. + * + * If epsilon is NULL, it means this is the first time Emg is being run, and we + * return g, initial viscosity. + */ + + /*output: */ + double viscosity3d; + + /*input strain rate: */ + double exx,eyy,exy,exz,eyz,ezz; + + /*Intermediary value A and exponent e: */ + double A,e; + double eps0; + + eps0=10^-27; + + if (n==1){ + /*Viscous behaviour! viscosity3d=B: */ + viscosity3d=B; + } + else{ + if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && + (epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){ + viscosity3d=pow(10,14); + } + else{ + + /*Retrive strain rate components: */ + exx=*(epsilon+0); + eyy=*(epsilon+1); + ezz=*(epsilon+2); //not used + exy=*(epsilon+3); + exz=*(epsilon+4); + eyz=*(epsilon+5); + + /*Build viscosity: viscosity3d=2*B/(2*A^e) */ + A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy+pow(eps0,2); + if(A==0){ + /*Maxiviscosity3dm viscosity for 0 shear areas: */ + viscosity3d=4.5*pow(10,17); + } + else{ + e=(n-1)/2/n; + + viscosity3d=2*B/(2*pow(A,e)); + } + } + } + #ifdef _DEBUG_ + _printf_("Viscosity %lf\n",viscosity3d); + #endif + + /*Assign output pointers:*/ + *pviscosity3d=viscosity3d; +} double Matice::GetB(){ return B; Index: /issm/trunk/src/c/objects/Matice.h =================================================================== --- /issm/trunk/src/c/objects/Matice.h (revision 393) +++ /issm/trunk/src/c/objects/Matice.h (revision 394) @@ -35,4 +35,5 @@ void GetViscosity2(double* pviscosity2, double* pepsilon); void GetViscosity3d(double* pviscosity3d, double* pepsilon); + void GetViscosity3dStokes(double* pviscosity3d, double* epsilon); Object* copy(); double GetB(); Index: /issm/trunk/src/c/objects/Pengrid.cpp =================================================================== --- /issm/trunk/src/c/objects/Pengrid.cpp (revision 393) +++ /issm/trunk/src/c/objects/Pengrid.cpp (revision 394) @@ -136,5 +136,5 @@ return my_rank; } -void Pengrid::DistributeNumDofs(int* numdofspernode,int analysis_type){return;} +void Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type){return;} #undef __FUNCT__ @@ -181,7 +181,57 @@ #define __FUNCT__ "Pengrid::PenaltyCreateKMatrix" void Pengrid::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){ - throw ErrorException(__FUNCT__," not implemented yet!"); -} - + + if ((analysis_type==DiagnosticStokesAnalysisEnum())){ + + PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type); + + } + else{ + throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet")); + } + +} + +#undef __FUNCT__ +#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixDiagnosticStokes" +void Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,void* vinputs,double kmax,int analysis_type){ + + const int numgrids=1; + const int NDOF3=3; + const int numdof=numgrids*NDOF3; + int doflist[numdof]; + int numberofdofspernode; + + int dofs1=0; + int dofs2=1; + double slope[2]; + int found=0; + double Ke[4][4]={{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}}; + + ParameterInputs* inputs=NULL; + + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /*Get dof list: */ + GetDofList(&doflist[0],&numberofdofspernode); + + /*recover slope: */ + found=inputs->Recover("bedslopex",&slope[0],1,&dofs1,numgrids,(void**)&node); + if(!found)throw ErrorException(__FUNCT__," bedslopex needed in inputs!"); + found=inputs->Recover("bedslopey",&slope[1],1,&dofs2,numgrids,(void**)&node); + if(!found)throw ErrorException(__FUNCT__," bedslopey needed in inputs!"); + + //Create elementary matrix: add penalty to contrain wb (wb=ub*db/dx+vb*db/dy) + Ke[0][0]=-slope[0]*kmax*pow(10.0,penalty_offset); + Ke[1][1]=-slope[1]*kmax*pow(10.0,penalty_offset); + Ke[2][2]=kmax*pow(10,penalty_offset); + + + /*Add Ke to global matrix Kgg: */ + MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES); +} + #undef __FUNCT__ #define __FUNCT__ "Pengrid::PenaltyCreatePVector" @@ -194,2 +244,17 @@ } + +void Pengrid::GetDofList(int* doflist,int* pnumberofdofspernode){ + + int j; + int doflist_per_node[MAXDOFSPERNODE]; + int numberofdofspernode; + + node->GetDofList(&doflist_per_node[0],&numberofdofspernode); + for(j=0;jGetRhoWater(); + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /*recover extra inputs from users, at current convergence iteration: */ + inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); + + /*Initialize Ke_temp to zero befor adding terms */ + for (i=0;i<27;i++){ + for (j=0;j<27;j++){ + Ke_temp[i][j]=0; + } + } + + /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore + get tria gaussian points as well as segment gaussian points. For tria gaussian + points, order of integration is 2, because we need to integrate the product tB*D*B' + which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian + points, same deal, which yields 3 gaussian points.*/ + + + area_order=5; + num_vert_gauss=5; + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); + + /* Start looping on the number of gaussian points: */ + for (igarea=0; igareaGetViscosity3dStokes(&viscosity,&epsilon[0]); + + /*Get B and Bprime matrices: */ + GetBStokes(&B[0][0],&xyz_list[0][0],gauss_coord); + GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss_coord); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],&gauss_coord[0]); + + /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant + * onto this scalar matrix, so that we win some computational time: */ + D_scalar=gauss_weight*Jdet; + for (i=0;i<6;i++){ + D[i][i]=D_scalar*viscosity; + } + for (i=6;i<8;i++){ + D[i][i]=-D_scalar*stokesreconditioning; + } + + + /* Do the triple product tB*D*Bprime: */ + MatrixMultiply(&B[0][0],8,27,1,&D[0][0],8,8,0,&tBD[0][0],0); + MatrixMultiply(&tBD[0][0],27,8,0,&B_prime[0][0],8,27,0,&Ke_gaussian[0][0],0); + + /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */ + for(i=0;i<27;i++){ + for(j=0;j<27;j++){ + Ke_temp[i][j]+=Ke_gaussian[i][j]; + } + } + } + } + + if((onbed==1) && (shelf==0)){ + + + /*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=&vxvyvz_list[0][0]; + friction->p=p; + friction->q=q; + friction->analysis_type=analysis_type; + + /*Compute alpha2_list: */ + FrictionGetAlpha2(&alpha2_list[0],friction); + + /*Erase friction object: */ + DeleteFriction(&friction); + + for(i=0;iGetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria); + + /*Get L matrix if viscous basal drag present: */ + GetLStokes(&LStokes[0][0], gauss_coord_tria); + GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss_coord_tria, gauss_coord); + + /*Compute strain rate: */ + GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord); + + /*Get viscosity at last iteration: */ + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + + /*Get normal vecyor to the bed */ + SurfaceNormal(&surface_normal[0],xyz_list_tria); + + bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result + bed_normal[1]=-surface_normal[1]; + bed_normal[2]=-surface_normal[2]; + + /*Calculate DL on gauss point */ + tria->GetParameterValue(&alpha2_gauss,&alpha2_list[0],gauss_coord_tria); + + DLStokes[0][0]=alpha2_gauss*gauss_weight*Jdet2d; + DLStokes[1][1]=alpha2_gauss*gauss_weight*Jdet2d; + DLStokes[2][2]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[0]*bed_normal[2]; + DLStokes[3][3]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[1]*bed_normal[2]; + DLStokes[4][4]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[0]*bed_normal[2]; + DLStokes[5][5]=-alpha2_gauss*gauss_weight*Jdet2d*bed_normal[1]*bed_normal[2]; + DLStokes[6][6]=-viscosity*gauss_weight*Jdet2d*bed_normal[0]; + DLStokes[7][7]=-viscosity*gauss_weight*Jdet2d*bed_normal[1]; + DLStokes[8][8]=-viscosity*gauss_weight*Jdet2d*bed_normal[2]; + DLStokes[9][8]=-viscosity*gauss_weight*Jdet2d*bed_normal[0]/2.0; + DLStokes[10][10]=-viscosity*gauss_weight*Jdet2d*bed_normal[1]/2.0; + DLStokes[11][11]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[0]; + DLStokes[12][12]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[1]; + DLStokes[13][13]=stokesreconditioning*gauss_weight*Jdet2d*bed_normal[2]; + + /* Do the triple product tL*D*L: */ + MatrixMultiply(&LStokes[0][0],14,numdof2d,1,&DLStokes[0][0],14,14,0,&tLDStokes[0][0],0); + MatrixMultiply(&tLDStokes[0][0],numdof2d,14,0,&LprimeStokes[0][0],14,numdof2d,0,&Ke_drag_gaussian[0][0],0); + + for(i=0;iGetRhoWater(); + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + + /*recover extra inputs from users, at current convergence iteration: */ + inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /*Initialize Ke_temp to zero befor adding terms */ + for (i=0;i<27;i++){ + Pe_temp[i]=0; + for (j=0;j<3;j++){ + Ke_temp[i][j]=0; + } + } + + + + /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore + get tria gaussian points as well as segment gaussian points. For tria gaussian + points, order of integration is 2, because we need to integrate the product tB*D*B' + which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian + points, same deal, which yields 3 gaussian points.*/ + + + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); + + /* Start looping on the number of gaussian points: */ + for (igarea=0; igareaGetViscosity3dStokes(&viscosity,&epsilon[0]); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); + + /* Get nodal functions */ + GetNodalFunctionsStokes(&l1l7[0], gauss_coord); + + /* Build gaussian vector */ + for(i=0;iGetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria); + + /* Get bed at gaussian point */ + GetParameterValue(&bed,&b[0],gauss_coord); + + /*Get L matrix : */ + tria->GetL(&L[0], &xyz_list[0][0], gauss_coord_tria,1); + + /*Get water_pressure at gaussian point */ + water_pressure=gravity*rho_water*bed; + + /*Get normal vecyor to the bed */ + SurfaceNormal(&surface_normal[0],xyz_list_tria); + + bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result + bed_normal[1]=-surface_normal[1]; + bed_normal[2]=-surface_normal[2]; + + for(i=0;i=3)]); %only keep vz if running with more than 3 dofs per node + soln(count).u_g=get(inputs,'velocity',[1 1 (m.parameters.numberofdofspernode>=3) (m.parameters.numberofdofspernode>=3)]); %only keep vz if running with more than 3 dofs per node soln(count).u_f=[]; Index: /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp =================================================================== --- /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp (revision 393) +++ /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp (revision 394) @@ -39,5 +39,4 @@ /*Create loads: */ CreateLoads(&loads,model,MODEL); - /*Create parameters: */ CreateParameters(¶meters,model,MODEL);