Index: ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 25488) +++ ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 25489) @@ -35,6 +35,44 @@ /*Assign output pointers:*/ *pvector=vector; } /*}}}*/ +void GetVectoronBaseFromInputsx(Vector** pvector,FemModel* femmodel,int name,int type){ /*{{{*/ + + int i, domaintype; + Vector* vector=NULL; + + switch(type){ + case ElementSIdEnum: + vector=new Vector(femmodel->elements->NumberOfElements()); + break; + case VertexPIdEnum: case VertexSIdEnum: + vector=new Vector(femmodel->vertices->NumberOfVertices()); + break; + case NodesEnum:case NodeSIdEnum: + vector=new Vector(femmodel->nodes->NumberOfNodes()); + break; + default: + _error_("vector type: " << EnumToStringx(type) << " not supported yet!"); + } + + /*Look up in elements*/ + for(i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + element->GetVectorFromInputs(vector,name,type); + break; + case Domain3DEnum: + if(!element->IsOnBase()) continue; + element->GetVectorFromInputs(vector,name,type); + break; + default: _error_("mesh "<Assemble(); + /*Assign output pointers:*/ + *pvector=vector; +} /*}}}*/ void GetVectorFromInputsx(Vector** pvector,FemModel* femmodel,int name,int type,IssmDouble time){/*{{{*/ int i; @@ -78,18 +116,35 @@ /*Assign output pointers:*/ *pvector=vector; }/*}}}*/ +void GetVectoronBaseFromInputsx(IssmDouble** pvector,FemModel* femmodel,int name, int type){/*{{{*/ + + /*output: */ + IssmDouble* vector=NULL; + + /*intermediary: */ + Vector* vec_vector=NULL; + + GetVectoronBaseFromInputsx(&vec_vector,femmodel,name,type); + vector=vec_vector->ToMPISerial(); + + /*Free ressources:*/ + delete vec_vector; + + /*Assign output pointers:*/ + *pvector=vector; +}/*}}}*/ void GetVectorFromInputsx(IssmDouble** pvector,int* pvector_size, FemModel* femmodel,int name){ /*{{{*/ int interpolation_type; - /*this one is special: we don't specify the type, but let the nature of the inputs dictace. + /*this one is special: we don't specify the type, but let the nature of the inputs dictace. * P0 -> ElementSIdEnum, P1 ->VertexSIdEnum: */ /*We go find the input of the first element, and query its interpolation type: */ Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(0)); - Input* input=element->GetInput(name); + Input* input=element->GetInput(name); if (!input) _error_("could not find input: " << name); - interpolation_type=input->GetInputInterpolationType(); + interpolation_type=input->GetInputInterpolationType(); if(interpolation_type==P0Enum){ *pvector_size=femmodel->elements->NumberOfElements(); GetVectorFromInputsx(pvector,femmodel,name, ElementSIdEnum); Index: ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h =================================================================== --- ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h (revision 25488) +++ ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h (revision 25489) @@ -1,5 +1,5 @@ /*!\file: GetVectorFromInputsx.h - */ + */ #ifndef _GETVECTORFROMINPUTSXX_H #define _GETVECTORFROMINPUTSXX_H @@ -7,9 +7,12 @@ #include "../../classes/classes.h" /* local prototypes: */ -void GetVectorFromInputsx( Vector** pvector,FemModel* femmodel,int name,int type); -void GetVectorFromInputsx( IssmDouble** pvector,FemModel* femmodel,int name,int type); -void GetVectorFromInputsx(IssmDouble** pvector,int* pvector_size, FemModel* femmodel,int name); -void GetVectorFromInputsx(Vector** pvector,FemModel* femmodel,int name,int type,IssmDouble time); +void GetVectorFromInputsx(Vector** pvector,FemModel* femmodel,int name,int type); +void GetVectoronBaseFromInputsx(Vector** pvector,FemModel* femmodel,int name,int type); +void GetVectorFromInputsx(Vector** pvector,FemModel* femmodel,int name,int type,IssmDouble time); +void GetVectorFromInputsx(IssmDouble** pvector,FemModel* femmodel,int name,int type); +void GetVectoronBaseFromInputsx(IssmDouble** pvector,FemModel* femmodel,int name,int type); +void GetVectorFromInputsx(IssmDouble** pvector,int* pvector_size, FemModel* femmodel,int name); + #endif /* _GETVECTORFROMINPUTSXX_H */ Index: ../trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 25488) +++ ../trunk-jpl/src/c/cores/hydrology_core.cpp (revision 25489) @@ -115,10 +115,10 @@ /*Proceed now to heads computations*/ solutionsequence_hydro_nonlinear(femmodel); /*If we have a sub-timestep we store the substep inputs in a transient input here*/ - femmodel->StackTransientInputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput); + femmodel->StackTransientInputonBasex(&substepinput[0],&transientinput[0],subtime,numaveragedinput); } /*averaging the stack*/ - femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,hydro_averaging); + femmodel->AverageTransientInputonBasex(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,hydro_averaging); /*And reseting to global time*/ femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum); femmodel->parameters->SetParam(global_time,TimeEnum); Index: ../trunk-jpl/src/c/classes/FemModel.h =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.h (revision 25488) +++ ../trunk-jpl/src/c/classes/FemModel.h (revision 25489) @@ -178,7 +178,9 @@ void UpdateConstraintsL2ProjectionEPLx(IssmDouble* pL2count); void InitTransientInputx(int* transientinput_enum,int numoutputs); void StackTransientInputx(int* input_enum,int* transientinput_enum,IssmDouble hydrotime,int numoutputs); + void StackTransientInputonBasex(int* input_enum,int* transientinput_enum,IssmDouble hydrotime,int numoutputs); void AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs,int averaging_method); + void AverageTransientInputonBasex(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs,int averaging_method); void UpdateConstraintsx(void); int UpdateVertexPositionsx(void); Index: ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 25488) +++ ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 25489) @@ -687,8 +687,8 @@ bool element_active; Element* element=NULL; - - for(int i=0;ielements->Size();i++){ + int elementssize=femmodel->elements->Size(); + for(int i=0;i(femmodel->elements->GetObjectByOffset(i)); Input* input=element->GetInput(HydrologydcMaskEplactiveNodeEnum); _assert_(input); @@ -747,8 +747,8 @@ bool element_active; Element* element=NULL; - - for(int i=0;ielements->Size();i++){ + int elementssize = femmodel->elements->Size(); + for(int i=0;i(femmodel->elements->GetObjectByOffset(i)); Input* input=element->GetInput(HydrologydcMaskThawedNodeEnum); _assert_(input); Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 25488) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 25489) @@ -48,7 +48,6 @@ #include "../modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h" #include "../modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h" #include "../modules/NodalValuex/NodalValuex.h" -#include "../modules/GetVectorFromInputsx/GetVectorFromInputsx.h" #include "../modules/AverageOntoPartitionx/AverageOntoPartitionx.h" /*}}}*/ @@ -4542,9 +4541,9 @@ int npart; /*retrieve partition vectors for responses that are scaled:*/ - this->parameters->FindParam(&response_partitions,&response_partitions_num,NULL,NULL,QmuResponsePartitionsEnum); - this->parameters->FindParam(&response_partitions_npart,NULL,NULL,QmuResponsePartitionsNpartEnum); - + this->parameters->FindParam(&response_partitions,&response_partitions_num,NULL,NULL,QmuResponsePartitionsEnum); + this->parameters->FindParam(&response_partitions_npart,NULL,NULL,QmuResponsePartitionsNpartEnum); + /*retrieve my_rank: */ my_rank=IssmComm::GetRank(); @@ -5020,7 +5019,7 @@ recurence=new Vector(numnodes); this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); - GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum); + GetVectoronBaseFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum); for (int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); @@ -5067,13 +5066,25 @@ /*Update node activation accordingly*/ int counter = 0; //this is probably not acurate but we are only interested in positivity + Element* basalelement=NULL; for(int i=0;iSize();i++){ Element *element = xDynamicCast(elements->GetObjectByOffset(i)); - int numnodes = element->GetNumberOfNodes(); + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + break; + case Domain3DEnum: + if(!element->IsOnBase()) continue; + basalelement = element->SpawnBasalElement(); + break; + default: _error_("mesh "<GetNumberOfNodes(); IssmDouble *base = xNew(numnodes); - element->GetInputListOnNodes(&base[0],BaseEnum); + basalelement->GetInputListOnNodes(&base[0],BaseEnum); for(int in=0;inGetNode(in); + Node* node=basalelement->GetNode(in); if(serial_active[node->Sid()]==1.){ node->Activate(); if(!node->IsClone()) counter++; @@ -5084,6 +5095,10 @@ } } xDelete(base); + if(domaintype!=Domain2DhorizontalEnum){ + basalelement->DeleteMaterials(); + delete basalelement; + } } xDelete(serial_active); delete effanalysis; @@ -5133,7 +5148,7 @@ } /*for other cases we just grab the mask from the initialisation value*/ else{ - GetVectorFromInputsx(&serial_mask,this,HydrologydcMaskThawedNodeEnum,NodeSIdEnum); + GetVectoronBaseFromInputsx(&serial_mask,this,HydrologydcMaskThawedNodeEnum,NodeSIdEnum); } /*Update Mask and elementize*/ InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskThawedNodeEnum,NodeSIdEnum); @@ -5280,6 +5295,61 @@ } } /*}}}*/ +void FemModel::StackTransientInputonBasex(int* input_enum,int* transientinput_enum,IssmDouble subtime,int numoutputs){ /*{{{*/ + + Element* basalelement=NULL; + int domaintype; + this->parameters->FindParam(&domaintype,DomainTypeEnum); + + for(int i=0;iSize();j++){ + /*Get the right transient input*/ + Element* element=xDynamicCast(elements->GetObjectByOffset(j)); + /*Get basal element*/ + switch(domaintype){ + case Domain2DhorizontalEnum: + break; + case Domain3DEnum: + if(!element->IsOnBase()) continue; + break; + default: _error_("mesh "<inputs->GetTransientInput(transientinput_enum[i]); + const int numvertices = element->GetNumberOfVertices(); + IssmDouble* values = xNew(numvertices); + int *vertexlids = xNew(numvertices); + switch(domaintype){ + case Domain2DhorizontalEnum: + element->GetInputListOnVertices(&values[0],input_enum[i]); //this is the enum to stack + element->GetVerticesLidList(vertexlids); + transientinput->AddTriaTimeInput(subtime,numvertices,vertexlids,values,P1Enum); + break; + case Domain3DEnum:{ + element->GetInputListOnVertices(&values[0],input_enum[i]); //this is the enum to stack + element->GetVerticesLidList(vertexlids); + Penta* penta=xDynamicCast(elements->GetObjectByOffset(j)); + for(;;){ + transientinput->AddPentaTimeInput(subtime,numvertices,vertexlids,values,P1Enum); + if (penta->IsOnSurface()) break; + penta=penta->GetUpperPenta(); _assert_(penta->Id()!=element->id); + penta->GetVerticesLidList(vertexlids); + } + break; + } + default: _error_("mesh "<(values); + xDelete(vertexlids); + } + } + } +} +/*}}}*/ void FemModel::AverageTransientInputx(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs, int averaging_method){ /*{{{*/ for(int i=0;iparameters->FindParam(&domaintype,DomainTypeEnum); + + for(int i=0;ielements->Size();j++){ + Element* element = xDynamicCast(elements->GetObjectByOffset(j)); + /*Get basal element*/ + switch(domaintype){ + case Domain2DhorizontalEnum: + element->CreateInputTimeAverage(transientinput_enum[i],averagedinput_enum[i],init_time,end_time,averaging_method);; + break; + case Domain3DEnum: + if(!element->IsOnBase()) continue; + element->CreateInputTimeAverage(transientinput_enum[i],averagedinput_enum[i],init_time,end_time,averaging_method); + break; + default: _error_("mesh "<