Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5923) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5924) @@ -3002,12 +3002,5 @@ void Penta::CreatePVectorBalancedthickness( Vec pg){ - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*If on water, skip: */ - if(IsOnWater())return; - - /*Is this element on the bed? :*/ - if(!IsOnBed())return; + if (!IsOnBed() || IsOnWater()) return; /*Depth Averaging Vx and Vy*/ @@ -3015,8 +3008,10 @@ this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); - /*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); + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + ElementVector* pe=tria->CreatePVectorBalancedthickness(); delete tria->matice; delete tria; + pe->AddToGlobal(pg,NULL); + delete pe; /*Delete Vx and Vy averaged*/ @@ -3024,4 +3019,5 @@ this->inputs->DeleteInput(VyAverageEnum); + /*Clean up and return*/ return; } @@ -3030,12 +3026,5 @@ void Penta::CreatePVectorBalancedvelocities( Vec pg){ - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*If on water, skip: */ - if(IsOnWater())return; - - /*Is this element on the bed? :*/ - if(!IsOnBed())return; + if (!IsOnBed() || IsOnWater()) return; /*Depth Averaging Vx and Vy*/ @@ -3043,8 +3032,10 @@ this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); - /*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); + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + ElementVector* pe=tria->CreatePVectorBalancedvelocities(); delete tria->matice; delete tria; + pe->AddToGlobal(pg,NULL); + delete pe; /*Delete Vx and Vy averaged*/ @@ -3052,4 +3043,5 @@ this->inputs->DeleteInput(VyAverageEnum); + /*Clean up and return*/ return; } @@ -3729,12 +3721,5 @@ void Penta::CreatePVectorPrognostic( Vec pg){ - /*Collapsed formulation: */ - Tria* tria=NULL; - - /*If on water, skip: */ - if(IsOnWater())return; - - /*Is this element on the bed? :*/ - if(!IsOnBed())return; + if (!IsOnBed() || IsOnWater()) return; /*Depth Averaging Vx and Vy*/ @@ -3742,8 +3727,10 @@ this->InputDepthAverageAtBase(VyEnum,VyAverageEnum); - /*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); + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + ElementVector* pe=tria->CreatePVectorPrognostic(); delete tria->matice; delete tria; + pe->AddToGlobal(pg,NULL); + delete pe; /*Delete Vx and Vy averaged*/ @@ -3751,4 +3738,5 @@ this->inputs->DeleteInput(VyAverageEnum); + /*Clean up and return*/ return; } Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5923) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5924) @@ -752,33 +752,25 @@ case DiagnosticHorizAnalysisEnum: pe=CreatePVectorDiagnosticMacAyeal(); - pe->AddToGlobal(pg,pf); - delete pe; break; case AdjointHorizAnalysisEnum: pe=CreatePVectorAdjointHoriz(); - pe->AddToGlobal(pg,pf); - delete pe; break; case DiagnosticHutterAnalysisEnum: pe=CreatePVectorDiagnosticHutter(); - pe->AddToGlobal(pg,pf); - delete pe; break; case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: pe=CreatePVectorSlope(); - pe->AddToGlobal(pg,pf); - delete pe; break; case PrognosticAnalysisEnum: - CreatePVectorPrognostic(pg); + pe=CreatePVectorPrognostic(); break; case BalancedthicknessAnalysisEnum: - CreatePVectorBalancedthickness(pg); + pe=CreatePVectorBalancedthickness(); break; case AdjointBalancedthicknessAnalysisEnum: - CreatePVectorAdjointBalancedthickness(pg); + pe=CreatePVectorAdjointBalancedthickness(); break; case BalancedvelocitiesAnalysisEnum: - CreatePVectorBalancedvelocities( pg); + pe=CreatePVectorBalancedvelocities(); break; default: @@ -786,4 +778,9 @@ } + /*Add to global Vector*/ + if(pe){ + pe->AddToGlobal(pg,pf); + delete pe; + } } /*}}}*/ @@ -3451,13 +3448,12 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorBalancedthickness{{{1*/ -void Tria::CreatePVectorBalancedthickness(Vec pg){ +ElementVector* Tria::CreatePVectorBalancedthickness(void){ switch(GetElementType()){ case P1Enum: - CreatePVectorBalancedthickness_CG( pg); + return CreatePVectorBalancedthickness_CG(); break; case P1DGEnum: - CreatePVectorBalancedthickness_DG( pg); - break; + return CreatePVectorBalancedthickness_DG(); default: ISSMERROR("Element type %s not supported yet",EnumToString(GetElementType())); @@ -3466,5 +3462,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorBalancedthickness_CG{{{1*/ -void Tria::CreatePVectorBalancedthickness_CG(Vec pg ){ +ElementVector* Tria::CreatePVectorBalancedthickness_CG(void){ /*Constants*/ @@ -3477,12 +3473,12 @@ double dhdt_g,melting_g,accumulation_g,Jdettria; double L[NUMVERTICES]; - double pe_g[NUMVERTICES] ={0.0}; 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 inputs :*/ Input* accumulation_input=inputs->GetInput(AccumulationRateEnum); ISSMASSERT(accumulation_input); Input* melting_input=inputs->GetInput(MeltingRateEnum); ISSMASSERT(melting_input); @@ -3502,16 +3498,14 @@ GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); - for( i=0; iweight*(accumulation_g-melting_g-dhdt_g)*L[i]; - } - - VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); + for(i=0;ivalues[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g-dhdt_g)*L[i]; + } /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return pe; } /*}}}*/ /*FUNCTION Tria::CreatePVectorBalancedthickness_DG {{{1*/ -void Tria::CreatePVectorBalancedthickness_DG(Vec pg){ +ElementVector* Tria::CreatePVectorBalancedthickness_DG(void){ /*Constants*/ @@ -3524,12 +3518,12 @@ double melting_g,accumulation_g,dhdt_g,Jdettria; double L[NUMVERTICES]; - double pe_g[NUMVERTICES] ={0.0}; 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: */ Input* accumulation_input=inputs->GetInput(AccumulationRateEnum); ISSMASSERT(accumulation_input); Input* melting_input=inputs->GetInput(MeltingRateEnum); ISSMASSERT(melting_input); @@ -3549,16 +3543,14 @@ GetL(&L[0], &xyz_list[0][0], gauss,NDOF1); - for(i=0;iweight*(accumulation_g-melting_g-dhdt_g)*L[i]; - } - - VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); + for(i=0;ivalues[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g-dhdt_g)*L[i]; + } /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return pe; } /*}}}*/ /*FUNCTION Tria::CreatePVectorBalancedvelocities {{{1*/ -void Tria::CreatePVectorBalancedvelocities(Vec pg){ +ElementVector* Tria::CreatePVectorBalancedvelocities(void){ /*Constants*/ @@ -3571,11 +3563,12 @@ double Jdettria,accumulation_g,melting_g; double L[NUMVERTICES]; - double pe_g[NUMVERTICES]={0.0}; GaussTria* gauss=NULL; + /*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: */ Input* accumulation_input=inputs->GetInput(AccumulationRateEnum); ISSMASSERT(accumulation_input); Input* melting_input=inputs->GetInput(MeltingRateEnum); ISSMASSERT(melting_input); @@ -3593,13 +3586,10 @@ melting_input->GetParameterValue(&melting_g,gauss); - for( i=0; iweight*(accumulation_g-melting_g)*L[i]; - - } - - VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); + for(i=0;ivalues[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g)*L[i]; + } /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return pe; } /*}}}*/ @@ -3742,5 +3732,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorAdjointBalancedthickness{{{1*/ -void Tria::CreatePVectorAdjointBalancedthickness(Vec p_g){ +ElementVector* Tria::CreatePVectorAdjointBalancedthickness(void){ /*Constants*/ @@ -3755,12 +3745,12 @@ double l1l2l3[3]; double pe_g_gaussian[numdof]; - double pe_g[numdof] ={0.0}; 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 and parameters: */ Input* thickness_input =inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); Input* thicknessobs_input=inputs->GetInput(ThicknessObsEnum);ISSMASSERT(thicknessobs_input); @@ -3780,17 +3770,10 @@ weights_input->GetParameterValue(&weight, gauss); - for (i=0;iweight*l1l2l3[i]; - } - - /*Add pe_g_gaussian vector to pe_g: */ - for(i=0; ivalues[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*l1l2l3[i]; + } /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return pe; } /*}}}*/ @@ -4256,13 +4239,11 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorPrognostic{{{1*/ -void Tria::CreatePVectorPrognostic(Vec pg){ +ElementVector* Tria::CreatePVectorPrognostic(void){ switch(GetElementType()){ case P1Enum: - CreatePVectorPrognostic_CG( pg); - break; + return CreatePVectorPrognostic_CG(); case P1DGEnum: - CreatePVectorPrognostic_DG( pg); - break; + return CreatePVectorPrognostic_DG(); default: ISSMERROR("Element type %s not supported yet",EnumToString(GetElementType())); @@ -4271,5 +4252,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorPrognostic_CG {{{1*/ -void Tria::CreatePVectorPrognostic_CG(Vec pg){ +ElementVector* Tria::CreatePVectorPrognostic_CG(void){ /* local declarations */ @@ -4281,15 +4262,8 @@ int* doflist=NULL; int numberofdofspernode=1; - - /* gaussian points: */ int ig; GaussTria* gauss=NULL; - - /* matrix */ - double pe_g[NUMVERTICES]={0.0}; double L[NUMVERTICES]; double Jdettria; - - /*input parameters for structural analysis (diagnostic): */ double accumulation_g; double melting_g; @@ -4297,12 +4271,11 @@ double dt; - /*retrieve some parameters: */ + /*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); this->parameters->FindParam(&dt,DtEnum); - - /* 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* accumulation_input=inputs->GetInput(AccumulationRateEnum); ISSMASSERT(accumulation_input); Input* melting_input=inputs->GetInput(MeltingRateEnum); ISSMASSERT(melting_input); @@ -4315,51 +4288,32 @@ gauss->GaussPoint(ig); - /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); - - /*Get L matrix: */ GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); - /* Get accumulation, melting and thickness at gauss point */ accumulation_input->GetParameterValue(&accumulation_g,gauss); melting_input->GetParameterValue(&melting_g,gauss); thickness_input->GetParameterValue(&thickness_g,gauss); - /* Add value into pe_g: */ - for( i=0; iweight*(thickness_g+dt*(accumulation_g-melting_g))*L[i]; - - } - - /*Add pe_g to global matrix Kgg: */ - VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); + for(i=0;ivalues[i]+=Jdettria*gauss->weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i]; + } /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); - + return pe; } /*}}}*/ /*FUNCTION Tria::CreatePVectorPrognostic_DG {{{1*/ -void Tria::CreatePVectorPrognostic_DG(Vec pg){ +ElementVector* Tria::CreatePVectorPrognostic_DG(void){ /* local declarations */ + const int numdof=NDOF1*NUMVERTICES; int i,j; - - /* node data: */ - const int numdof=NDOF1*NUMVERTICES; + int ig; double xyz_list[NUMVERTICES][3]; int* doflist=NULL; int numberofdofspernode=1; - - /* gaussian points: */ - int ig; GaussTria* gauss=NULL; - - /* matrix */ - double pe_g[NUMVERTICES]={0.0}; double L[NUMVERTICES]; double Jdettria; - - /*input parameters for structural analysis (diagnostic): */ double accumulation_g; double melting_g; @@ -4367,12 +4321,11 @@ double dt; - /*retrieve some parameters: */ + /*Initialize Element vector and return if necessary*/ + if(IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ this->parameters->FindParam(&dt,DtEnum); - - /* 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* accumulation_input=inputs->GetInput(AccumulationRateEnum); ISSMASSERT(accumulation_input); Input* melting_input=inputs->GetInput(MeltingRateEnum); ISSMASSERT(melting_input); @@ -4385,26 +4338,17 @@ gauss->GaussPoint(ig); - /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss); - - /*Get L matrix: */ GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode); - /* Get accumulation, melting and thickness at gauss point */ accumulation_input->GetParameterValue(&accumulation_g,gauss); melting_input->GetParameterValue(&melting_g,gauss); thickness_input->GetParameterValue(&thickness_g,gauss); - /* Add value into pe_g: */ - for( i=0; iweight*(thickness_g+dt*(accumulation_g-melting_g))*L[i]; - - } - - /*Add pe_g to global matrix Kgg: */ - VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES); + for(i=0;ivalues[i]+=Jdettria*gauss->weight*(thickness_g+dt*(accumulation_g-melting_g))*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 5923) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5924) @@ -138,17 +138,17 @@ ElementMatrix* CreateKMatrixSlope(void); ElementMatrix* CreateKMatrixThermal(void); - void CreatePVectorBalancedthickness(Vec pg); - void CreatePVectorBalancedthickness_DG(Vec pg); - void CreatePVectorBalancedthickness_CG(Vec pg); - void CreatePVectorBalancedvelocities(Vec pg); + ElementVector* CreatePVectorBalancedthickness(void); + ElementVector* CreatePVectorBalancedthickness_DG(void); + ElementVector* CreatePVectorBalancedthickness_CG(void); + ElementVector* CreatePVectorBalancedvelocities(void); void CreatePVectorDiagnosticBaseVert(Vec pg); ElementVector* CreatePVectorDiagnosticMacAyeal(void); ElementVector* CreatePVectorAdjointHoriz(void); void CreatePVectorAdjointStokes(Vec pg); - void CreatePVectorAdjointBalancedthickness(Vec pg); + ElementVector* CreatePVectorAdjointBalancedthickness(void); ElementVector* CreatePVectorDiagnosticHutter(void); - void CreatePVectorPrognostic(Vec pg); - void CreatePVectorPrognostic_CG(Vec pg); - void CreatePVectorPrognostic_DG(Vec pg); + ElementVector* CreatePVectorPrognostic(void); + ElementVector* CreatePVectorPrognostic_CG(void); + ElementVector* CreatePVectorPrognostic_DG(void); ElementVector* CreatePVectorSlope(void); void CreatePVectorThermalSheet( Vec pg);