Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 1046) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 1047) @@ -4,5 +4,5 @@ #ifdef HAVE_CONFIG_H -#include "config.h" + #include "config.h" #else #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" @@ -33,9 +33,9 @@ Tria::Tria(int tria_id,int tria_mid,int tria_mparid,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,double tria_meanvel,double tria_epsvel, - double tria_viscosity_overshoot,int tria_artdiff){ - + double tria_accumulation[3],double tria_geothermalflux[3],int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel, + double tria_viscosity_overshoot,int tria_artdiff){ + int i; - + id=tria_id; mid=tria_mid; @@ -66,5 +66,5 @@ viscosity_overshoot=tria_viscosity_overshoot; artdiff=tria_artdiff; - + return; } @@ -156,8 +156,8 @@ /*get enum type of Tria: */ enum_type=TriaEnum(); - + /*marshall enum: */ memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); - + /*marshall Tria data: */ memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id); @@ -187,37 +187,37 @@ memcpy(marshalled_dataset,&viscosity_overshoot,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot); memcpy(marshalled_dataset,&artdiff,sizeof(artdiff));marshalled_dataset+=sizeof(artdiff); - + *pmarshalled_dataset=marshalled_dataset; return; } - + int Tria::MarshallSize(){ return sizeof(id) - +sizeof(mid) - +sizeof(mparid) - +sizeof(node_ids) - +sizeof(nodes) - +sizeof(node_offsets) - +sizeof(matice) - +sizeof(matice_offset) - +sizeof(matpar) - +sizeof(matpar_offset) - +sizeof(h) - +sizeof(s) - +sizeof(b) - +sizeof(k) - +sizeof(melting) - +sizeof(accumulation) - +sizeof(geothermalflux) - +sizeof(friction_type) - +sizeof(onbed) - +sizeof(p) - +sizeof(q) - +sizeof(shelf) - +sizeof(meanvel) - +sizeof(epsvel) - +sizeof(viscosity_overshoot) - +sizeof(artdiff) - +sizeof(int); //sizeof(int) for enum type + +sizeof(mid) + +sizeof(mparid) + +sizeof(node_ids) + +sizeof(nodes) + +sizeof(node_offsets) + +sizeof(matice) + +sizeof(matice_offset) + +sizeof(matpar) + +sizeof(matpar_offset) + +sizeof(h) + +sizeof(s) + +sizeof(b) + +sizeof(k) + +sizeof(melting) + +sizeof(accumulation) + +sizeof(geothermalflux) + +sizeof(friction_type) + +sizeof(onbed) + +sizeof(p) + +sizeof(q) + +sizeof(shelf) + +sizeof(meanvel) + +sizeof(epsvel) + +sizeof(viscosity_overshoot) + +sizeof(artdiff) + +sizeof(int); //sizeof(int) for enum type } @@ -291,5 +291,5 @@ int i; - + DataSet* loadsin=NULL; DataSet* nodesin=NULL; @@ -303,5 +303,5 @@ /*Link this element with its nodes, ie find pointers to the nodes in the nodes dataset.: */ ResolvePointers((Object**)nodes,node_ids,node_offsets,3,nodesin); - + /*Same for materials: */ ResolvePointers((Object**)&matice,&mid,&matice_offset,1,materialsin); @@ -317,9 +317,9 @@ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ if (analysis_type==ControlAnalysisEnum()){ - + CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type); } else if (analysis_type==DiagnosticAnalysisEnum()){ - + if (sub_analysis_type==HorizAnalysisEnum()){ @@ -361,5 +361,5 @@ int doflist[numdof]; int numberofdofspernode; - + /* gaussian points: */ int num_gauss,ig; @@ -375,5 +375,5 @@ double newviscosity; //viscosity double oldviscosity; //viscosity - + /* strain rate: */ double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ @@ -389,7 +389,7 @@ double Ke_gg[numdof][numdof]; //local element stiffness matrix double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. - + double Jdet; - + /*input parameters for structural analysis (diagnostic): */ double vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; @@ -414,10 +414,10 @@ for(i=0;iGetRhoIce()); printf(" gravity: %g \n",matpar->GetG()) - printf(" rho_water: %g \n",matpar->GetRhoWater()); + printf(" rho_water: %g \n",matpar->GetRhoWater()); printf(" Velocity: \n"); for (i=0;iGetViscosity2d(&viscosity, &epsilon[0]); matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]); - + /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); /* 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: */ + onto this scalar matrix, so that we win some computational time: */ newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); D_scalar=newviscosity*thickness*gauss_weight*Jdet; @@ -476,6 +476,6 @@ D[i][i]=D_scalar; } - -#ifdef _DEBUGELEMENTS_ + + #ifdef _DEBUGELEMENTS_ if(my_rank==RANK && id==ELID){ printf(" gaussian loop %i\n",ig); @@ -492,5 +492,5 @@ printf(" gauss_weight: %g \n",gauss_weight); } -#endif + #endif /*Get B and Bprime matrices: */ @@ -500,12 +500,12 @@ /* Do the triple product tB*D*Bprime: */ TripleMultiply( &B[0][0],3,numdof,1, - &D[0][0],3,3,0, - &Bprime[0][0],3,numdof,0, - &Ke_gg_gaussian[0][0],0); + &D[0][0],3,3,0, + &Bprime[0][0],3,numdof,0, + &Ke_gg_gaussian[0][0],0); /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ for( i=0; iRecover("dt",&dt); if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); @@ -645,5 +645,5 @@ v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]); v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]); - + K[0][0]=pow(Jdettria,.5)/2.0*fabs(v_gauss[0]); K[1][1]=pow(Jdettria,.5)/2.0*fabs(v_gauss[1]); @@ -671,16 +671,16 @@ /* Do the triple product tL*D*L: */ TripleMultiply( &L[0],1,numdof,1, - &DL_scalar,1,1,0, - &L[0],1,numdof,0, - &Ke_gg_gaussian[0][0],0); - + &DL_scalar,1,1,0, + &L[0],1,numdof,0, + &Ke_gg_gaussian[0][0],0); + /*Get B and B prime matrix: */ GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3); GetBPrime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); - + //Get vx, vy and their derivatives at gauss point GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3); GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3); - + GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3); GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3); @@ -702,12 +702,12 @@ TripleMultiply( &B[0][0],2,numdof,1, - &DL[0][0],2,2,0, - &B[0][0],2,numdof,0, - &Ke_gg_thickness1[0][0],0); + &DL[0][0],2,2,0, + &B[0][0],2,numdof,0, + &Ke_gg_thickness1[0][0],0); TripleMultiply( &B[0][0],2,numdof,1, - &DLprime[0][0],2,2,0, - &Bprime[0][0],2,numdof,0, - &Ke_gg_thickness2[0][0],0); + &DLprime[0][0],2,2,0, + &Bprime[0][0],2,numdof,0, + &Ke_gg_thickness2[0][0],0); /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ @@ -715,7 +715,7 @@ for( i=0; ielement_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char)); strcpy(friction->element_type,"2d"); - + friction->gravity=matpar->GetG(); friction->rho_ice=matpar->GetRhoIce(); @@ -974,14 +974,14 @@ DeleteFriction(&friction); -#ifdef _DEBUGELEMENTS_ + #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 + #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_ + #ifdef _DEBUGELEMENTS_ if(my_rank==RANK && id==ELID){ printf(" gaussian points: \n"); @@ -990,5 +990,5 @@ } } -#endif + #endif /* Start looping on the number of gaussian points: */ @@ -1025,5 +1025,5 @@ DL[i][i]=DL_scalar; } - + /* Do the triple producte tL*D*L: */ TripleMultiply( &L[0][0],2,numdof,1, @@ -1039,5 +1039,5 @@ MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); -cleanup_and_return: + cleanup_and_return: xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); @@ -1062,5 +1062,5 @@ int doflist[numdof]; int numberofdofspernode; - + /* gaussian points: */ int num_gauss,ig; @@ -1079,7 +1079,7 @@ double Ke_gg[numdof][numdof]; //local element stiffness matrix double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. - + double Jdet; - + /* Get node coordinates and dof list: */ GetElementNodeData( &xyz_list[0][0], nodes, numgrids); @@ -1100,5 +1100,5 @@ gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); - + /*Get L matrix: */ GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1); @@ -1106,5 +1106,5 @@ /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - + DL_scalar=gauss_weight*Jdet; @@ -1121,6 +1121,6 @@ /*Add Ke_gg to global matrix Kgg: */ MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES); - -cleanup_and_return: + + cleanup_and_return: xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); @@ -1132,21 +1132,21 @@ #define __FUNCT__ "Tria::CreatePVector" void Tria::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){ - + /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ if (analysis_type==ControlAnalysisEnum()){ - + CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); - + } else if (analysis_type==DiagnosticAnalysisEnum()){ if (sub_analysis_type==HorizAnalysisEnum()){ - + CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); - + } else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet")); } else if (analysis_type==SlopeComputeAnalysisEnum()){ - + CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type); } @@ -1174,5 +1174,5 @@ int doflist[numdof]; int numberofdofspernode; - + /* parameters: */ double plastic_stress; @@ -1215,5 +1215,5 @@ -#ifdef _DEBUGELEMENTS_ + #ifdef _DEBUGELEMENTS_ if(my_rank==RANK && id==ELID){ printf("gravity %g\n",matpar->GetG()); @@ -1224,5 +1224,5 @@ printf("drag [%g,%g,%g]\n",k[0],k[1],k[2]); } -#endif + #endif @@ -1230,5 +1230,5 @@ 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_ + #ifdef _DEBUGELEMENTS_ if(my_rank==RANK && id==ELID){ printf(" gaussian points: \n"); @@ -1237,5 +1237,5 @@ } } -#endif + #endif @@ -1251,7 +1251,7 @@ /*Compute thickness at gaussian point: */ GetParameterValue(&thickness, &h[0],gauss_l1l2l3); - + GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3); - + /*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the * element itself: */ @@ -1262,6 +1262,6 @@ /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /*Get nodal functions: */ + + /*Get nodal functions: */ GetNodalFunctions(l1l2l3, gauss_l1l2l3); @@ -1270,5 +1270,5 @@ -#ifdef _DEBUGELEMENTS_ + #ifdef _DEBUGELEMENTS_ if(my_rank==RANK && id==ELID){ printf(" gaussian %i\n",ig); @@ -1280,5 +1280,5 @@ if(friction_type==1)printf(" plastic_stress(%g)\n",plastic_stress); } -#endif + #endif /*Build pe_g_gaussian vector: */ @@ -1303,5 +1303,5 @@ } //for (ig=0; igterms\n",ig); @@ -1315,10 +1315,10 @@ } } -#endif + #endif /*Add pe_g to global vector pg: */ VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); -cleanup_and_return: + cleanup_and_return: xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); @@ -1342,5 +1342,5 @@ int doflist[numdof]; int numberofdofspernode; - + /* gaussian points: */ int num_gauss,ig; @@ -1391,9 +1391,9 @@ GetParameterDerivativeValue(&slope[0], ¶m[0],&xyz_list[0][0], gauss_l1l2l3); - + /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /*Get nodal functions: */ + + /*Get nodal functions: */ GetNodalFunctions(l1l2l3, gauss_l1l2l3); @@ -1414,5 +1414,5 @@ VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); -cleanup_and_return: + cleanup_and_return: xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); @@ -1445,5 +1445,5 @@ inputs->Recover("accumulation",&accumulation[0],1,dofs,3,(void**)nodes); inputs->Recover("geothermalflux",&geothermalflux[0],1,dofs,3,(void**)nodes); - + //Update material if necessary if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,3,(void**)nodes)){ @@ -1452,5 +1452,5 @@ matice->SetB(B_average); } - + if(inputs->Recover("B",&B_list[0],1,dofs,3,(void**)nodes)){ B_average=(B_list[0]+B_list[1]+B_list[2])/3.0; @@ -1459,5 +1459,5 @@ } - + void Tria::GetDofList(int* doflist,int* pnumberofdofspernode){ @@ -1465,5 +1465,5 @@ int doflist_per_node[MAXDOFSPERNODE]; int numberofdofspernode; - + for(i=0;i<3;i++){ nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode); @@ -1491,8 +1491,8 @@ #define __FUNCT__ "Tria::GetParameterValue" void Tria::GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3){ - + /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian * point specifie by gauss_l1l2l3: */ - + /*nodal functions: */ double l1l2l3[3]; @@ -1513,5 +1513,5 @@ #define __FUNCT__ "Tria::GetParameterDerivativeValue" void Tria::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3){ - + const int NDOF2=2; const int numgrids=3; @@ -1523,5 +1523,5 @@ * p is a vector of size 2x1 already allocated. */ - + double dh1dh2dh3_basic[NDOF2][numgrids]; //nodal derivative functions in basic coordinate system. @@ -1548,5 +1548,5 @@ GetB(&B[0][0], xyz_list, gauss_l1l2l3); -#ifdef _DEBUG_ + #ifdef _DEBUG_ printf("B for grid1 : [ %lf %lf \n",B[0][0],B[0][1]); printf(" [ %lf %lf \n",B[1][0],B[1][1]); @@ -1560,14 +1560,14 @@ printf(" [ %lf %lf \n",B[1][4],B[1][5]); printf(" [ %lf %lf ]\n",B[2][4],B[2][5]); - + for (i=0;i