Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 2710) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 2711) @@ -16,11 +16,14 @@ #include "../include/typedefs.h" +/*FUNCTION Penta constructor {{{1*/ Penta::Penta(){ return; } +/*}}}*/ +/*FUNCTION Penta creation {{{1*/ Penta::Penta( int penta_id, int penta_mid, int penta_mparid, int penta_numparid, int penta_node_ids[6], double penta_h[6], double penta_s[6], double penta_b[6], double penta_k[6], int penta_friction_type, - double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6], - int penta_thermal_steadystate,bool penta_onwater){ - + double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6], + int penta_thermal_steadystate,bool penta_onwater){ + int i; @@ -61,8 +64,11 @@ return; } - +/*}}}*/ +/*FUNCTION Penta destructor {{{1*/ Penta::~Penta(){ return; } +/*}}}*/ +/*FUNCTION Penta Echo {{{1*/ void Penta::Echo(void){ @@ -82,5 +88,5 @@ printf(" b=[%g,%g,%g,%g,%g,%g]\n",b[0],b[1],b[2],b[3],b[4],b[5]); printf(" k=[%g,%g,%g,%g,%g,%g]\n",k[0],k[1],k[2],k[3],k[4],k[5]); - + printf(" friction_type: %i\n",friction_type); printf(" p: %g\n",p); @@ -91,5 +97,5 @@ printf(" onsurface: %i\n",onsurface); printf(" collapse: %i\n",collapse); - + printf(" melting=[%g,%g,%g,%g,%g,%g]\n",melting[0],melting[1],melting[2],melting[3],melting[4],melting[5]); printf(" accumulation=[%g,%g,%g,%g,%g,%g]\n",accumulation[0],accumulation[1],accumulation[2],accumulation[3],accumulation[4],accumulation[5]); @@ -98,4 +104,6 @@ return; } +/*}}}*/ +/*FUNCTION Penta DeepEcho {{{1*/ void Penta::DeepEcho(void){ @@ -115,5 +123,5 @@ printf(" b=[%i,%i,%i,%i,%i,%i]\n",b[0],b[1],b[2],b[3],b[4],b[5]); printf(" k=[%i,%i,%i,%i,%i,%i]\n",k[0],k[1],k[2],k[3],k[4],k[5]); - + printf(" friction_type: %i\n",friction_type); printf(" p: %g\n",p); @@ -124,5 +132,5 @@ printf(" onsurface: %i\n",onsurface); printf(" collapse: %i\n",collapse); - + printf(" melting=[%i,%i,%i,%i,%i,%i]\n",melting[0],melting[1],melting[2],melting[3],melting[4],melting[5]); printf(" accumulation=[%i,%i,%i,%i,%i,%i]\n",accumulation[0],accumulation[1],accumulation[2],accumulation[3],accumulation[4],accumulation[5]); @@ -131,4 +139,6 @@ return; } +/*}}}*/ +/*FUNCTION Penta Marshall {{{1*/ void Penta::Marshall(char** pmarshalled_dataset){ @@ -141,8 +151,8 @@ /*get enum type of Penta: */ enum_type=PentaEnum(); - + /*marshall enum: */ memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); - + /*marshall Penta data: */ memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id); @@ -175,47 +185,50 @@ memcpy(marshalled_dataset,&geothermalflux,sizeof(geothermalflux));marshalled_dataset+=sizeof(geothermalflux); memcpy(marshalled_dataset,&thermal_steadystate,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate); - + *pmarshalled_dataset=marshalled_dataset; return; } - +/*}}}*/ +/*FUNCTION Penta MarshallSize {{{1*/ int Penta::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(numparid)+ - +sizeof(numpar)+ - +sizeof(numpar_offset)+ - sizeof(h)+ - sizeof(s)+ - sizeof(b)+ - sizeof(k)+ - sizeof(friction_type)+ - sizeof(p)+ - sizeof(q)+ - sizeof(shelf)+ - sizeof(onbed)+ - sizeof(onwater)+ - sizeof(onsurface)+ - sizeof(collapse)+ - sizeof(melting)+ - sizeof(accumulation)+ - sizeof(geothermalflux)+ - sizeof(thermal_steadystate) + - 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(numparid)+ + +sizeof(numpar)+ + +sizeof(numpar_offset)+ + sizeof(h)+ + sizeof(s)+ + sizeof(b)+ + sizeof(k)+ + sizeof(friction_type)+ + sizeof(p)+ + sizeof(q)+ + sizeof(shelf)+ + sizeof(onbed)+ + sizeof(onwater)+ + sizeof(onsurface)+ + sizeof(collapse)+ + sizeof(melting)+ + sizeof(accumulation)+ + sizeof(geothermalflux)+ + sizeof(thermal_steadystate) + + sizeof(int); //sizeof(int) for enum type +} +/*}}}*/ +/*FUNCTION Penta GetName {{{1*/ char* Penta::GetName(void){ return "penta"; } - +/*}}}*/ +/*FUNCTION Penta Demarshall {{{1*/ void Penta::Demarshall(char** pmarshalled_dataset){ @@ -269,4 +282,6 @@ return; } +/*}}}*/ +/*FUNCTION Penta Enum {{{1*/ int Penta::Enum(void){ @@ -274,18 +289,23 @@ } -int Penta::GetId(void){ return id; } - +/*}}}*/ +/*FUNCTION Penta GetId {{{1*/ +int Penta::GetId(void){ + return id; +} +/*}}}*/ +/*FUNCTION Penta MyRank {{{1*/ int Penta::MyRank(void){ extern int my_rank; return my_rank; } - - -#undef __FUNCT__ +/*}}}*/ +/*FUNCTION Penta Configure {{{1*/ +#undef __FUConfigure NCT__ #define __FUNCT__ "Penta::Configure" void Penta::Configure(void* ploadsin,void* pnodesin,void* pmaterialsin,void* pparametersin){ int i; - + DataSet* loadsin=NULL; DataSet* nodesin=NULL; @@ -301,5 +321,5 @@ /*Link this element with its nodes, ie find pointers to the nodes in the nodes dataset.: */ ResolvePointers((Object**)nodes,node_ids,node_offsets,6,nodesin); - + /*Same for materials: */ ResolvePointers((Object**)&matice,&mid,&matice_offset,1,materialsin); @@ -310,5 +330,6 @@ } - +/*}}}*/ +/*FUNCTION Penta CreateKMatrix {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::CreateKMatrix" @@ -323,34 +344,34 @@ } else if (analysis_type==DiagnosticAnalysisEnum()){ - + if (sub_analysis_type==HorizAnalysisEnum()){ - + CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type); } else if (sub_analysis_type==VertAnalysisEnum()){ - + CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type,sub_analysis_type); } else if (sub_analysis_type==StokesAnalysisEnum()){ - + CreateKMatrixDiagnosticStokes( Kgg,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()){ - + CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type); } else if (analysis_type==PrognosticAnalysisEnum()){ - + CreateKMatrixPrognostic( Kgg,inputs,analysis_type,sub_analysis_type); } else if (analysis_type==ThermalAnalysisEnum()){ - + CreateKMatrixThermal( Kgg,inputs,analysis_type,sub_analysis_type); } else if (analysis_type==MeltingAnalysisEnum()){ - + CreateKMatrixMelting( Kgg,inputs,analysis_type,sub_analysis_type); } @@ -360,5 +381,6 @@ } - +/*}}}*/ +/*FUNCTION Penta CreateKMatrixDiagnosticHoriz {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta:CreateKMatrixDiagnosticHoriz" @@ -375,6 +397,6 @@ int doflist[numdof]; int numberofdofspernode; - - + + /* 3d gaussian points: */ int num_gauss,ig; @@ -392,5 +414,5 @@ int num_area_gauss; double gauss_weight; - + /* 2d gaussian point: */ int num_gauss2d; @@ -424,5 +446,5 @@ double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag double Jdet; - + /*slope: */ double slope[2]={0.0,0.0}; @@ -491,5 +513,5 @@ } - #ifdef _DEBUGELEMENTS_ +#ifdef _DEBUGELEMENTS_ if(my_rank==RANK && id==ELID){ printf("El id %i Rank %i PentaElement input list before gaussian loop: \n",ELID,RANK); @@ -508,5 +530,5 @@ printf(" temperature_average [%g %g %g %g %g %g]\n",temperature_average_list[0],temperature_average_list[1],temperature_average_list[2],temperature_average_list[3],temperature_average_list[4],temperature_average_list[5]); } - #endif +#endif /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore @@ -521,5 +543,5 @@ GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); - #ifdef _DEBUGGAUSS_ +#ifdef _DEBUGGAUSS_ if(my_rank==RANK && id==ELID){ printf("El id %i Rank %i PentaElement gauss points\n",ELID,RANK); @@ -531,15 +553,15 @@ } } - #endif +#endif /* Start looping on the number of gaussian points: */ for (ig1=0; ig1GetViscosity3d(&viscosity, &epsilon[0]); @@ -562,8 +584,8 @@ /* Get Jacobian determinant: */ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4); - + /*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: */ - + newviscosity=viscosity+numpar->viscosity_overshoot*(viscosity-oldviscosity); D_scalar=newviscosity*gauss_weight*Jdet; @@ -571,10 +593,10 @@ D[i][i]=D_scalar; } - + /* Do the triple product tB*D*Bprime: */ TripleMultiply( &B[0][0],5,numdof,1, - &D[0][0],5,5,0, - &Bprime[0][0],5,numdof,0, - &Ke_gg_gaussian[0][0],0); + &D[0][0],5,5,0, + &Bprime[0][0],5,numdof,0, + &Ke_gg_gaussian[0][0],0); /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ @@ -586,9 +608,9 @@ } //for (ig2=0; ig2GetRhoWater(); @@ -861,8 +885,8 @@ /* 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 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; @@ -888,5 +912,5 @@ /*Compute strain rate: */ GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord); - + /*Get viscosity: */ matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); @@ -926,9 +950,9 @@ /*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(); @@ -953,5 +977,5 @@ } } - + xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); xfree((void**)&third_gauss_area_coord); xfree((void**)&area_gauss_weights); GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2); @@ -964,25 +988,25 @@ gauss_coord[2]=*(third_gauss_area_coord+igarea); gauss_coord[3]=-1; - + gauss_coord_tria[0]=*(first_gauss_area_coord+igarea); gauss_coord_tria[1]=*(second_gauss_area_coord+igarea); gauss_coord_tria[2]=*(third_gauss_area_coord+igarea); - + /*Get the Jacobian determinant */ tria->GetJacobianDeterminant3d(&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]; @@ -1018,5 +1042,5 @@ } } //if ( (onbed==1) && (shelf==0)) - + /*Reduce the matrix */ ReduceMatrixStokes(&Ke_reduced[0][0], &Ke_temp[0][0]); @@ -1030,7 +1054,7 @@ /*Add Ke_gg to global matrix Kgg: */ MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES); - - cleanup_and_return: - + +cleanup_and_return: + /*Free ressources:*/ xfree((void**)&first_gauss_area_coord); @@ -1043,5 +1067,6 @@ return; } - +/*}}}*/ +/*FUNCTION Penta CreatePVector {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::CreatePVector" @@ -1050,9 +1075,9 @@ /*Just branch to the correct element stiffness matrix 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()){ @@ -1060,9 +1085,9 @@ } else if (sub_analysis_type==VertAnalysisEnum()){ - + CreatePVectorDiagnosticVert( pg,inputs,analysis_type,sub_analysis_type); } else if (sub_analysis_type==StokesAnalysisEnum()){ - + CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type); } @@ -1070,17 +1095,17 @@ } else if (analysis_type==SlopeComputeAnalysisEnum()){ - + CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type); } else if (analysis_type==PrognosticAnalysisEnum()){ - + CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type); } else if (analysis_type==ThermalAnalysisEnum()){ - + CreatePVectorThermal( pg,inputs,analysis_type,sub_analysis_type); } else if (analysis_type==MeltingAnalysisEnum()){ - + CreatePVectorMelting( pg,inputs,analysis_type,sub_analysis_type); } @@ -1090,4 +1115,6 @@ } +/*}}}*/ +/*FUNCTION Penta UpdateFromInputs {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::UpdateFromInputs" @@ -1126,5 +1153,5 @@ } } - + if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,6,(void**)nodes)){ if(matice && collapse){ @@ -1136,5 +1163,5 @@ } } - + if(inputs->Recover("B",&B_list[0],1,dofs,6,(void**)nodes)){ if(matice){ @@ -1145,14 +1172,16 @@ } - +/*}}}*/ +/*FUNCTION Penta GetMatPar {{{1*/ void* Penta::GetMatPar(){ return matpar; } - +/*}}}*/ +/*FUNCTION Penta GetShelf {{{1*/ int Penta::GetShelf(){ return shelf; } - - +/*}}}*/ +/*FUNCTION Penta GetNodes {{{1*/ void Penta::GetNodes(void** vpnodes){ int i; @@ -1163,25 +1192,31 @@ } } - +/*}}}*/ +/*FUNCTION Penta GetOnBed {{{1*/ int Penta::GetOnBed(){ return onbed; } - -void Penta::GetThicknessList(double* thickness_list){ +/*}}}*/ +/*FUNCTION Penta GetThicknessList {{{1*/ +void Penta::GetThicknessList(double* thickness_list){ int i; for(i=0;i<6;i++)thickness_list[i]=h[i]; } -void Penta::GetBedList(double* bed_list){ - +/*}}}*/ +/*FUNCTION Penta GetBedList {{{1*/ +void Penta::GetBedList(double* bed_list){ + int i; for(i=0;i<6;i++)bed_list[i]=b[i]; } - +/*}}}*/ +/*FUNCTION Penta copy {{{1*/ Object* Penta::copy() { return new Penta(*this); } - +/*}}}*/ +/*FUNCTION Penta Du {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::Du" @@ -1193,5 +1228,5 @@ /*If on water, skip: */ if(onwater)return; - + /*Bail out if this element if: * -> Non collapsed and not on the surface @@ -1210,5 +1245,5 @@ } else{ - + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). tria->Du(du_g,inputs,analysis_type,sub_analysis_type); @@ -1217,9 +1252,10 @@ } } - +/*}}}*/ +/*FUNCTION Penta Gradj {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::Gradj" void Penta::Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){ - + /*If on water, skip grad (=0): */ if(onwater)return; @@ -1233,9 +1269,10 @@ else throw ErrorException(__FUNCT__,exprintf("%s%s","control type not supported yet: ",control_type)); } - +/*}}}*/ +/*FUNCTION Penta GradjDrag {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::GradjDrag" void Penta::GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type){ - + Tria* tria=NULL; @@ -1267,5 +1304,6 @@ else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet")); } - +/*}}}*/ +/*FUNCTION Penta GradjB {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::GradjB" @@ -1296,12 +1334,13 @@ } } - +/*}}}*/ +/*FUNCTION Penta Misfit {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::Misfit" double Penta::Misfit(void* inputs,int analysis_type,int sub_analysis_type){ - + double J; Tria* tria=NULL; - + /*If on water, return 0: */ if(onwater)return 0; @@ -1330,5 +1369,6 @@ } } - +/*}}}*/ +/*FUNCTION Penta SpawnTria {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::SpawnTria" @@ -1345,5 +1385,5 @@ double tria_accumulation[3]; double tria_geothermalflux[3]; - + /*configuration: */ int tria_node_ids[3]; @@ -1370,5 +1410,5 @@ tria_melting[1]=melting[g1]; tria_melting[2]=melting[g2]; - + tria_accumulation[0]=accumulation[g0]; tria_accumulation[1]=accumulation[g1]; @@ -1401,6 +1441,6 @@ } - - +/*}}}*/ +/*FUNCTION Penta GetDofList {{{1*/ void Penta::GetDofList(int* doflist,int* pnumberofdofspernode){ @@ -1408,5 +1448,5 @@ int doflist_per_node[MAXDOFSPERNODE]; int numberofdofspernode; - + for(i=0;i<6;i++){ nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode); @@ -1420,5 +1460,6 @@ } - +/*}}}*/ +/*FUNCTION Penta GetDofList1 {{{1*/ void Penta::GetDofList1(int* doflist){ @@ -1429,4 +1470,6 @@ } +/*}}}*/ +/*FUNCTION Penta GetStrainRate {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::GetStrainRate" @@ -1442,5 +1485,5 @@ GetB(&B[0][0], xyz_list, gauss_l1l2l3l4); - #ifdef _ISSM_DEBUG_ +#ifdef _ISSM_DEBUG_ printf("B for grid1 : [ %lf %lf \n",B[0][0],B[0][1]); printf(" [ %lf %lf \n",B[1][0],B[1][1]); @@ -1448,5 +1491,5 @@ printf(" [ %lf %lf ]\n",B[3][0],B[3][1]); printf(" [ %lf %lf ]\n",B[4][0],B[4][1]); - + printf("B for grid2 : [ %lf %lf \n",B[0][2],B[0][3]); printf(" [ %lf %lf \n",B[1][2],B[1][3]); @@ -1454,5 +1497,5 @@ printf(" [ %lf %lf ]\n",B[3][2],B[3][3]); printf(" [ %lf %lf ]\n",B[4][2],B[4][3]); - + printf("B for grid3 : [ %lf %lf \n", B[0][4],B[0][5]); printf(" [ %lf %lf \n", B[1][4],B[1][5]); @@ -1460,5 +1503,5 @@ printf(" [ %lf %lf ]\n",B[3][4],B[3][5]); printf(" [ %lf %lf ]\n",B[4][4],B[4][5]); - + printf("B for grid4 : [ %lf %lf \n", B[0][6],B[0][7]); printf(" [ %lf %lf \n", B[1][6],B[1][7]); @@ -1466,5 +1509,5 @@ printf(" [ %lf %lf ]\n",B[3][6],B[3][7]); printf(" [ %lf %lf ]\n",B[4][6],B[4][7]); - + printf("B for grid5 : [ %lf %lf \n", B[0][8],B[0][9]); printf(" [ %lf %lf \n", B[1][8],B[1][9]); @@ -1479,13 +1522,14 @@ printf(" [ %lf %lf ]\n",B[4][10],B[4][11]); - #endif +#endif /*Multiply B by velocity, to get strain rate: */ MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0, - velocity,NDOF2*numgrids,1,0, - epsilon,0); - -} - + velocity,NDOF2*numgrids,1,0, + epsilon,0); + +} +/*}}}*/ +/*FUNCTION Penta GetB {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::GetB" @@ -1504,5 +1548,5 @@ * We assume B has been allocated already, of size: 5x(NDOF2*numgrids) */ - + int i; const int numgrids=6; @@ -1515,9 +1559,9 @@ GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4); - #ifdef _ISSM_DEBUG_ +#ifdef _ISSM_DEBUG_ for (i=0;iRecover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity!"); - inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes); + inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes); /*Get gaussian points and weights :*/ @@ -2266,5 +2319,5 @@ GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); - #ifdef _ISSM_DEBUG_ +#ifdef _ISSM_DEBUG_ for (i=0;iGetViscosity3dStokes(&viscosity,&epsilon[0]); - + /* Get Jacobian determinant: */ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); @@ -3203,5 +3263,5 @@ GetBStokes(&B[0][0],&xyz_list[0][0],gauss_coord); GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss_coord); - + /*Get bubble part of Bprime */ for(i=0;i<8;i++){ @@ -3242,5 +3302,5 @@ } } - + xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); xfree((void**)&third_gauss_area_coord); xfree((void**)&area_gauss_weights); GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2); @@ -3253,5 +3313,5 @@ gauss_coord[2]=*(third_gauss_area_coord+igarea); gauss_coord[3]=-1; - + gauss_coord_tria[0]=*(first_gauss_area_coord+igarea); gauss_coord_tria[1]=*(second_gauss_area_coord+igarea); @@ -3260,5 +3320,5 @@ /*Get the Jacobian determinant */ tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria); - + /* Get bed at gaussian point */ GetParameterValue(&bed,&b[0],gauss_coord); @@ -3266,5 +3326,5 @@ /*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; @@ -3272,5 +3332,5 @@ /*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]; @@ -3284,5 +3344,5 @@ } } //if ( (onbed==1) && (shelf==1)) - + /*Reduce the matrix */ ReduceVectorStokes(&Pe_reduced[0], &Ke_temp[0][0], &Pe_temp[0]); @@ -3304,9 +3364,10 @@ } - +/*}}}*/ +/*FUNCTION Penta ReduceVectorStokes {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::ReduceVectorStokes" void Penta::ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp){ - + int i,j; @@ -3347,9 +3408,10 @@ } } - +/*}}}*/ +/*FUNCTION Penta GetNodalFunctionsStokes {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::GetNodalFunctionsStokes" void Penta::GetNodalFunctionsStokes(double* l1l7, double* gauss_coord){ - + /*This routine returns the values of the nodal functions at the gaussian point.*/ @@ -3376,9 +3438,10 @@ } - +/*}}}*/ +/*FUNCTION Penta CreateKMatrixThermal {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::CreateKMatrixThermal" void Penta::CreateKMatrixThermal(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ - + /* local declarations */ int i,j; @@ -3466,5 +3529,5 @@ GetElementNodeData( &xyz_list[0][0], nodes, numgrids); GetDofList(&doflist[0],&numberofdofspernode); - + // /*recovre material parameters: */ rho_water=matpar->GetRhoWater(); @@ -3473,5 +3536,5 @@ heatcapacity=matpar->GetHeatCapacity(); thermalconductivity=matpar->GetThermalConductivity(); - + /*recover extra inputs from users, dt and velocity: */ found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); @@ -3488,9 +3551,9 @@ /* 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 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: */ area_order=2; @@ -3601,5 +3664,5 @@ } } - + /*Add Ke_gaussian to pKe: */ for(i=0;iCreateKMatrixMelting(Kgg,inputs, analysis_type,sub_analysis_type); @@ -3787,5 +3855,6 @@ } } - +/*}}}*/ +/*FUNCTION Penta CreatePVectorThermal {{{1*/ #undef __FUNCT__ #define __FUNCT__ "Penta::CreatePVectorThermal" @@ -3805,5 +3874,5 @@ /*Grid data: */ double xyz_list[numgrids][3]; - + /* gaussian points: */ int num_area_gauss,igarea,igvert; @@ -3818,5 +3887,5 @@ int area_order=2; int num_vert_gauss=3; - + double dt; double vx_list[numgrids]; @@ -3834,5 +3903,5 @@ /* element parameters: */ int friction_type; - + int dofs[3]={0,1,2}; int dofs1[1]={0}; @@ -3872,5 +3941,5 @@ GetElementNodeData( &xyz_list[0][0], nodes, numgrids); GetDofList(&doflist[0],&numberofdofspernode); - + /*recovre material parameters: */ rho_water=matpar->GetRhoWater(); @@ -3891,5 +3960,5 @@ if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!"); } - + for(i=0;i