Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5921) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5922) @@ -3771,17 +3771,14 @@ void Penta::CreatePVectorSlope( Vec pg){ - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*If on water, skip: */ - if(IsOnWater())return; - - /*Is this element on the bed? :*/ - if(!IsOnBed())return; - - /*Spawn Tria element from the base of the Penta: */ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg,NULL); + if (!IsOnBed() || IsOnWater()) return; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + ElementVector* pe=tria->CreatePVectorSlope(); delete tria->matice; delete tria; + pe->AddToGlobal(pg,NULL); + delete pe; + + /*clean up and return*/ return; } Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5921) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5922) @@ -761,8 +761,12 @@ break; case DiagnosticHutterAnalysisEnum: - CreatePVectorDiagnosticHutter( pg); + pe=CreatePVectorDiagnosticHutter(); + pe->AddToGlobal(pg,pf); + delete pe; break; case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: - CreatePVectorSlope( pg); + pe=CreatePVectorSlope(); + pe->AddToGlobal(pg,pf); + delete pe; break; case PrognosticAnalysisEnum: @@ -3687,5 +3691,5 @@ /*Initialize Element vector and return if necessary*/ if(IsOnWater()) return NULL; - ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters); + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum); /*Retrieve all inputs and parameters*/ @@ -4197,5 +4201,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorDiagnosticHutter{{{1*/ -void Tria::CreatePVectorDiagnosticHutter(Vec pg){ +ElementVector* Tria::CreatePVectorDiagnosticHutter(void){ /*Collapsed formulation: */ @@ -4205,5 +4209,4 @@ double constant_part,ub,vb; double rho_ice,gravity,n,B; - double pe_g[numdofs]; double slope[2]; double thickness; @@ -4212,16 +4215,13 @@ GaussTria* gauss=NULL; - /*If on water, skip: */ - if(IsOnWater())return; - - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /* Get parameters */ + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ rho_ice=matpar->GetRhoIce(); gravity=matpar->GetG(); n=matice->GetN(); B=matice->GetBbar(); - - /* Get slopes and thickness */ Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); ISSMASSERT(slopex_input); Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); ISSMASSERT(slopey_input); @@ -4241,21 +4241,16 @@ slope2=pow(slope[0],2)+pow(slope[1],2); - //compute constant_part constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2)); - //compute ub ub=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[0]; vb=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[1]; - pe_g[2*i]=(ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/(double)connectivity; - pe_g[2*i+1]=(vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/(double)connectivity; - - } - - VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES); + pe->values[2*i] =(ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/(double)connectivity; + pe->values[2*i+1]=(vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/(double)connectivity; + } /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return pe; } /*}}}*/ @@ -4415,40 +4410,27 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorSlope {{{1*/ -void Tria::CreatePVectorSlope( Vec pg){ - - int i,j; +ElementVector* Tria::CreatePVectorSlope(void){ /* node data: */ const int numdof=NDOF1*NUMVERTICES; + int i,j; + int ig; double xyz_list[NUMVERTICES][3]; int* doflist=NULL; - - /* gaussian points: */ - int ig; GaussTria* gauss=NULL; - - /* Jacobian: */ double Jdet; - - /*nodal functions: */ double l1l2l3[3]; - - /*element vector at the gaussian points: */ - double pe_g[numdof]={0.0}; double pe_g_gaussian[numdof]; double slope[2]; int analysis_type; - /*inputs :*/ + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); Input* slope_input=NULL; - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /*Retrieve all inputs we will be needing: */ if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==SurfaceSlopeYAnalysisEnum)){ slope_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(slope_input); @@ -4481,14 +4463,11 @@ /*Add pe_g_gaussian vector to pe_g: */ - for( i=0; ivalues[i]+=pe_g_gaussian[i]; + + } /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return pe; } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 5921) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5922) @@ -147,9 +147,9 @@ void CreatePVectorAdjointStokes(Vec pg); void CreatePVectorAdjointBalancedthickness(Vec pg); - void CreatePVectorDiagnosticHutter(Vec pg); + ElementVector* CreatePVectorDiagnosticHutter(void); void CreatePVectorPrognostic(Vec pg); void CreatePVectorPrognostic_CG(Vec pg); void CreatePVectorPrognostic_DG(Vec pg); - void CreatePVectorSlope( Vec pg); + ElementVector* CreatePVectorSlope(void); void CreatePVectorThermalSheet( Vec pg); void CreatePVectorThermalShelf( Vec pg);