Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5932) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5933) @@ -752,49 +752,31 @@ case DiagnosticHorizAnalysisEnum: pe=CreatePVectorDiagnosticHoriz(); - if(pe) pe->AddToGlobal(pg,pf); - if(pe) delete pe; break; case AdjointHorizAnalysisEnum: pe=CreatePVectorAdjointHoriz(); - if(pe) pe->AddToGlobal(pg,pf); - if(pe) delete pe; break; case DiagnosticHutterAnalysisEnum: pe=CreatePVectorDiagnosticHutter(); - if(pe) pe->AddToGlobal(pg,pf); - if(pe) delete pe; break; case DiagnosticVertAnalysisEnum: - CreatePVectorDiagnosticVert( pg); + pe=CreatePVectorDiagnosticVert(); break; case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: pe=CreatePVectorSlope(); - if(pe) pe->AddToGlobal(pg,pf); - if(pe) delete pe; break; case PrognosticAnalysisEnum: pe=CreatePVectorPrognostic(); - if(pe) pe->AddToGlobal(pg,pf); - if(pe) delete pe; break; case BalancedthicknessAnalysisEnum: pe=CreatePVectorBalancedthickness(); - if(pe) pe->AddToGlobal(pg,pf); - if(pe) delete pe; break; case BalancedvelocitiesAnalysisEnum: pe=CreatePVectorBalancedvelocities(); - if(pe) pe->AddToGlobal(pg,pf); - if(pe) delete pe; break; case ThermalAnalysisEnum: pe=CreatePVectorThermal(); - if(pe) pe->AddToGlobal(pg,pf); - if(pe) delete pe; break; case MeltingAnalysisEnum: pe=CreatePVectorMelting(); - if(pe) pe->AddToGlobal(pg,pf); - if(pe) delete pe; break; default: @@ -802,4 +784,9 @@ } + /*Add to global Vector*/ + if(pe){ + pe->AddToGlobal(pg,pf); + delete pe; + } } /*}}}*/ @@ -3058,5 +3045,5 @@ ElementVector* Penta::CreatePVectorCouplingPattynStokes(void){ - /*compute all load vectorsfor this element*/ + /*compute all load vectors for this element*/ ElementVector* pe1=CreatePVectorCouplingPattynStokesViscous(); ElementVector* pe2=CreatePVectorCouplingPattynStokesFriction(); @@ -3253,5 +3240,5 @@ ElementVector* Penta::CreatePVectorDiagnosticMacAyealPattyn(void){ - /*compute all load vectorsfor this element*/ + /*compute all load vectors for this element*/ ElementVector* pe1=CreatePVectorDiagnosticMacAyeal(); ElementVector* pe2=CreatePVectorDiagnosticPattyn(); @@ -3267,5 +3254,5 @@ ElementVector* Penta::CreatePVectorDiagnosticPattynStokes(void){ - /*compute all load vectorsfor this element*/ + /*compute all load vectors for this element*/ ElementVector* pe1=CreatePVectorDiagnosticPattyn(); ElementVector* pe2=CreatePVectorDiagnosticStokes(); @@ -3435,5 +3422,5 @@ ElementVector* Penta::CreatePVectorDiagnosticStokes(void){ - /*compute all load vectorsfor this element*/ + /*compute all load vectors for this element*/ ElementVector* pe1=CreatePVectorDiagnosticStokesViscous(); ElementVector* pe2=CreatePVectorDiagnosticStokesShelf(); @@ -3642,55 +3629,42 @@ /*}}}*/ /*FUNCTION Penta::CreatePVectorDiagnosticVert {{{1*/ -void Penta::CreatePVectorDiagnosticVert( Vec pg){ - - int i; - - /* node data: */ +ElementVector* Penta::CreatePVectorDiagnosticVert(void){ + + /*compute all load vectors for this element*/ + ElementVector* pe1=CreatePVectorDiagnosticVertVolume(); + ElementVector* pe2=CreatePVectorDiagnosticVertBase(); + ElementVector* pe =new ElementVector(pe1,pe2); + + /*clean-up and return*/ + delete pe1; + delete pe2; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticVertVolume {{{1*/ +ElementVector* Penta::CreatePVectorDiagnosticVertVolume(void){ + const int numdof=NDOF1*NUMVERTICES; double xyz_list[NUMVERTICES][3]; int* doflist=NULL; - - /* gaussian points: */ + int i; int ig; GaussPenta *gauss=NULL; - - /* Jacobian: */ double Jdet; - - /*element vector at the gaussian points: */ - double pe_g[numdof]={0.0}; - double pe_g_gaussian[numdof]; double l1l6[6]; - - /*Spawning: */ Tria* tria=NULL; - - /*input parameters for structural analysis (diagnostic): */ double du[3]; double dv[3]; double dw[3]; double dudx,dvdy,dwdz; - - /*Get some parameters*/ int approximation; - /*If on water, skip: */ - if(IsOnWater())return; + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); inputs->GetParameterValue(&approximation,ApproximationEnum); - - /*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the - *diagnostic base vertical stifness: */ - if(IsOnBed()){ - tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock - tria->CreatePVectorDiagnosticBaseVert(pg); - delete tria->matice; delete tria; - } - - /*Now, handle the standard penta element: */ - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /*Retrieve all inputs we will be needing: */ Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); @@ -3706,5 +3680,4 @@ gauss->GaussPoint(ig); - /*Get velocity derivative, with respect to x and y: */ vx_input->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss); vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss); @@ -3717,24 +3690,27 @@ dvdy=dv[1]; - /* Get Jacobian determinant: */ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); - - /*Get nodal functions: */ GetNodalFunctionsP1(l1l6, gauss); - /*Build pe_g_gaussian vector: */ - for (i=0;iweight*l1l6[i]; - } - - /*Add pe_g_gaussian vector to pe_g: */ - for( i=0; ivalues[i] += (dudx+dvdy+dwdz)*Jdet*gauss->weight*l1l6[i]; + } + + /*Clean up and return*/ delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticVertBase {{{1*/ +ElementVector* Penta::CreatePVectorDiagnosticVertBase(void){ + + if (!IsOnBed() || IsOnWater()) return NULL; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + ElementVector* pe=tria->CreatePVectorDiagnosticBaseVert(); + delete tria->matice; delete tria; + + /*Clean up and return*/ + return pe; } /*}}}*/ @@ -3783,5 +3759,5 @@ ElementVector* Penta::CreatePVectorThermal(void){ - /*compute all load vectorsfor this element*/ + /*compute all load vectors for this element*/ ElementVector* pe1=CreatePVectorThermalVolume(); ElementVector* pe2=CreatePVectorThermalSheet(); Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 5932) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 5933) @@ -168,5 +168,7 @@ ElementVector* CreatePVectorDiagnosticStokesShelf(void); ElementVector* CreatePVectorAdjointStokes(void); - void CreatePVectorDiagnosticVert( Vec pg); + ElementVector* CreatePVectorDiagnosticVert(void); + ElementVector* CreatePVectorDiagnosticVertVolume(void); + ElementVector* CreatePVectorDiagnosticVertBase(void); ElementVector* CreatePVectorMelting(void); ElementVector* CreatePVectorPrognostic(void); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5932) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5933) @@ -3595,5 +3595,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorDiagnosticBaseVert {{{1*/ -void Tria::CreatePVectorDiagnosticBaseVert(Vec pg){ +ElementVector* Tria::CreatePVectorDiagnosticBaseVert(void){ /*Constants*/ @@ -3603,5 +3603,4 @@ int i,j,ig; int approximation; - int* doflist=NULL; double xyz_list[NUMVERTICES][3]; double Jdet; @@ -3609,13 +3608,12 @@ double slope[2]; double L[NUMVERTICES]; - double pe_g[numdof]={0.0}; - double pe_g_gaussian[numdof]; GaussTria* gauss=NULL; - /* Get node coordinates and dof list: */ + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,NoneApproximationEnum,GsetEnum); - - /*Retrieve all inputs we will be needing: */ inputs->GetParameterValue(&approximation,ApproximationEnum); Input* bed_input=inputs->GetInput(BedEnum); ISSMASSERT(bed_input); @@ -3649,16 +3647,10 @@ GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); - for(i=0;iweight*(vx*dbdx+vy*dbdy-vz-meltingvalue)*L[i]; - } - - for( i=0; ivalues[i]+=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-meltingvalue)*L[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 5932) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5933) @@ -142,5 +142,5 @@ ElementVector* CreatePVectorBalancedthickness_CG(void); ElementVector* CreatePVectorBalancedvelocities(void); - void CreatePVectorDiagnosticBaseVert(Vec pg); + ElementVector* CreatePVectorDiagnosticBaseVert(void); ElementVector* CreatePVectorDiagnosticMacAyeal(void); ElementVector* CreatePVectorAdjointHoriz(void);