Index: /issm/trunk/src/c/objects/Loads/Icefront.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Icefront.cpp (revision 5934) +++ /issm/trunk/src/c/objects/Loads/Icefront.cpp (revision 5935) @@ -325,30 +325,24 @@ /*Retrieve parameters: */ + ElementVector* pe=NULL; int analysis_type; this->parameters->FindParam(&analysis_type,AnalysisTypeEnum); /*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */ - if (analysis_type==DiagnosticHorizAnalysisEnum){ - int type; - inputs->GetParameterValue(&type,TypeEnum); - if (type==MacAyeal2dIceFrontEnum){ - CreatePVectorDiagnosticMacAyeal2d( pg,pf); - } - else if (type==MacAyeal3dIceFrontEnum){ - CreatePVectorDiagnosticMacAyeal3d( pg); - } - else if (type==PattynIceFrontEnum){ - CreatePVectorDiagnosticPattyn( pg); - } - else if (type==StokesIceFrontEnum){ - CreatePVectorDiagnosticStokes( pg); - } - else ISSMERROR("Icefront type %s not supported yet",EnumToString(type)); - } - else if (analysis_type==AdjointHorizAnalysisEnum){ - return; - } - else{ - ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); + switch(analysis_type){ + case DiagnosticHorizAnalysisEnum: + pe=CreatePVectorDiagnosticHoriz(); + break; + case AdjointHorizAnalysisEnum: + pe=CreatePVectorAdjointHoriz(); + break; + default: + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); + } + + /*Add to global Vector*/ + if(pe){ + pe->AddToGlobal(pg,pf); + delete pe; } } @@ -424,13 +418,37 @@ /*Icefront numerics: */ +/*FUNCTION Icefront::CreatePVectorDiagnosticHoriz {{{1*/ +ElementVector* Icefront::CreatePVectorDiagnosticHoriz(void){ + + int type; + inputs->GetParameterValue(&type,TypeEnum); + + switch(type){ + case MacAyeal2dIceFrontEnum: + return CreatePVectorDiagnosticMacAyeal2d(); + case MacAyeal3dIceFrontEnum: + return CreatePVectorDiagnosticMacAyeal3d(); + case PattynIceFrontEnum: + return CreatePVectorDiagnosticPattyn(); + case StokesIceFrontEnum: + return CreatePVectorDiagnosticStokes(); + default: + ISSMERROR("Icefront type %s not supported yet",EnumToString(type)); + } +} +/*}}}*/ +/*FUNCTION Icefront::CreatePVectorAdjointHoriz {{{1*/ +ElementVector* Icefront::CreatePVectorAdjointHoriz(void){ + + /*No load vector applied to the adjoint*/ + return NULL; +} +/*}}}*/ /*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal2d{{{1*/ -void Icefront::CreatePVectorDiagnosticMacAyeal2d( Vec pg,Vec pf){ +ElementVector* Icefront::CreatePVectorDiagnosticMacAyeal2d(void){ /*Constants*/ const int numnodes= NUMVERTICESSEG; const int numdofs = numnodes *NDOF2; - - /*output: */ - ElementVector* pe=NULL; /*Intermediary*/ @@ -440,5 +458,4 @@ double water_pressure,air_pressure,surface_under_water,base_under_water; double xyz_list[numnodes][3]; - double pe_g[numdofs] = {0.0}; double normal[2]; double L[2]; @@ -447,18 +464,19 @@ Tria* tria=((Tria*)element); - /* Get dof list and node coordinates: */ + /*Initialize Element vector and return if necessary*/ + if(tria->IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICESSEG,this->parameters); + + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes); - - /*Recover inputs and parameters */ Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); Input* bed_input =tria->inputs->GetInput(BedEnum); ISSMASSERT(bed_input); inputs->GetParameterValue(&fill,FillEnum); - rho_water=matpar->GetRhoWater(); rho_ice =matpar->GetRhoIce(); gravity =matpar->GetG(); - GetSegmentNormal(&normal[0],xyz_list); + /*Start looping on Gaussian points*/ index1=tria->GetNodeIndex(nodes[0]); index2=tria->GetNodeIndex(nodes[1]); @@ -492,27 +510,16 @@ for (int i=0;iweight*normal[0]*L[i]; - pe_g[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*L[i]; + pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*L[i]; + pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*L[i]; } } - /*Initialize element matrix: */ - pe=NewElementVector(nodes,NUMVERTICESSEG,this->parameters,MacAyealApproximationEnum); - - /*Add pe_g values to pe element stifness load: */ - pe->AddValues(&pe_g[0]); - - /*Add pe element load vector to global load vector: */ - pe->AddToGlobal(pg,pf); - - /*Free ressources:*/ - delete pe; - - /*Free ressources:*/ + /*Clean up and return*/ delete gauss; + return pe; } /*}}}*/ /*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal3d{{{1*/ -void Icefront::CreatePVectorDiagnosticMacAyeal3d( Vec pg){ +ElementVector* Icefront::CreatePVectorDiagnosticMacAyeal3d(void){ Icefront* icefront=NULL; @@ -525,6 +532,5 @@ /*Return if not on bed*/ - penta->inputs->GetParameterValue(&onbed,ElementOnBedEnum); - if(!onbed) return; + if(!penta->IsOnBed() || penta->IsOnWater()) return NULL; /*Spawn Tria and call MacAyeal2d*/ @@ -533,14 +539,15 @@ icefront->element=tria; icefront->inputs->AddInput(new IntInput(TypeEnum,MacAyeal2dIceFrontEnum)); - icefront->CreatePVectorDiagnosticMacAyeal2d(pg,NULL); - - /*clean-up*/ + ElementVector* pe=icefront->CreatePVectorDiagnosticMacAyeal2d(); + + /*clean-up and return*/ delete tria->matice; delete tria; delete icefront; + return pe; } /*}}}*/ /*FUNCTION Icefront::CreatePVectorDiagnosticPattyn{{{1*/ -void Icefront::CreatePVectorDiagnosticPattyn( Vec pg){ +ElementVector* Icefront::CreatePVectorDiagnosticPattyn(void){ /*Constants*/ @@ -555,16 +562,15 @@ double xyz_list[NUMVERTICESQUA][3]; double normal[3]; - double pe_g[numdofs] = {0.0}; double l1l4[4]; - int *doflist = NULL; GaussPenta *gauss = NULL; - Penta *penta = NULL; - - /* Get dof list and node coordinates: */ - penta=(Penta*)element; - GetDofList(&doflist,PattynApproximationEnum,GsetEnum); + + Penta* penta=(Penta*)element; + + /*Initialize Element vector and return if necessary*/ + if(penta->IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICESQUA,this->parameters); + + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESQUA); - - /*Recover inputs and parameters */ Input* surface_input =penta->inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input); inputs->GetParameterValue(&fill,FillEnum); @@ -608,20 +614,14 @@ pressure = ice_pressure + water_pressure + air_pressure; - for(i=0;iweight*pressure*l1l4[i]*normal[j]; - } - } - } - - VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES); - - /*Free ressources:*/ + for(i=0;ivalues[i*NDOF2+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j]; + } + + /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return pe; } /*}}}*/ /*FUNCTION Icefront::CreatePVectorDiagnosticStokes{{{1*/ -void Icefront::CreatePVectorDiagnosticStokes( Vec pg){ +ElementVector* Icefront::CreatePVectorDiagnosticStokes(void){ /*Constants*/ @@ -636,16 +636,15 @@ double xyz_list[NUMVERTICESQUA][3]; double normal[3]; - double pe_g[numdofs] = {0.0}; double l1l4[4]; - int *doflist = NULL; GaussPenta *gauss = NULL; - Penta *penta = NULL; - - /* Get dof list and node coordinates: */ - penta=(Penta*)element; - GetDofList(&doflist,StokesApproximationEnum,GsetEnum); + + Penta* penta=(Penta*)element; + + /*Initialize Element vector and return if necessary*/ + if(penta->IsOnWater()) return NULL; + ElementVector* pe=NewElementVector(nodes,NUMVERTICESQUA,this->parameters); + + /*Retrieve all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESQUA); - - /*Recover inputs and parameters */ inputs->GetParameterValue(&fill,FillEnum); rho_water=matpar->GetRhoWater(); @@ -687,15 +686,13 @@ for(i=0;iweight*pressure*l1l4[i]*normal[j]; - else pe_g[i*NDOF4+j]+=0; //pressure term + if(j<3) pe->values[i*NDOF4+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j]; + else pe->values[i*NDOF4+j]+=0; //pressure term } } } - VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES); - - /*Free ressources:*/ + /*Clean up and return*/ delete gauss; - xfree((void**)&doflist); + return pe; } /*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Icefront.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Icefront.h (revision 5934) +++ /issm/trunk/src/c/objects/Loads/Icefront.h (revision 5935) @@ -77,8 +77,10 @@ /*}}}*/ /*Load management: {{{1*/ - void CreatePVectorDiagnosticMacAyeal2d(Vec pg,Vec pf); - void CreatePVectorDiagnosticMacAyeal3d(Vec pg); - void CreatePVectorDiagnosticPattyn(Vec pg); - void CreatePVectorDiagnosticStokes(Vec pg); + ElementVector* CreatePVectorAdjointHoriz(void); + ElementVector* CreatePVectorDiagnosticHoriz(void); + ElementVector* CreatePVectorDiagnosticMacAyeal2d(void); + ElementVector* CreatePVectorDiagnosticMacAyeal3d(void); + ElementVector* CreatePVectorDiagnosticPattyn(void); + ElementVector* CreatePVectorDiagnosticStokes(void); void GetDofList(int** pdoflist,int approximation_enum,int setenum); void GetSegmentNormal(double* normal,double xyz_list[2][3]);