Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 25488) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 25489) @@ -688,6 +688,6 @@ 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)); @@ -748,6 +748,6 @@ 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)); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 25488) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 25489) @@ -49,5 +49,4 @@ #include "../modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h" #include "../modules/NodalValuex/NodalValuex.h" -#include "../modules/GetVectorFromInputsx/GetVectorFromInputsx.h" #include "../modules/AverageOntoPartitionx/AverageOntoPartitionx.h" /*}}}*/ @@ -4543,7 +4542,7 @@ /*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(); @@ -5021,5 +5020,5 @@ 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++){ @@ -5068,11 +5067,23 @@ /*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(); @@ -5085,4 +5096,8 @@ } xDelete(base); + if(domaintype!=Domain2DhorizontalEnum){ + basalelement->DeleteMaterials(); + delete basalelement; + } } xDelete(serial_active); @@ -5134,5 +5149,5 @@ /*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*/ @@ -5281,4 +5296,59 @@ } /*}}}*/ +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){ /*{{{*/ @@ -5287,4 +5357,26 @@ Element* element = xDynamicCast(elements->GetObjectByOffset(j)); element->CreateInputTimeAverage(transientinput_enum[i],averagedinput_enum[i],init_time,end_time,averaging_method); + } + } +}/*}}}*/ +void FemModel::AverageTransientInputonBasex(int* transientinput_enum,int* averagedinput_enum,IssmDouble init_time,IssmDouble end_time,int numoutputs, int averaging_method){ /*{{{*/ + + int domaintype; + this->parameters->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 "<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); Index: /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 25488) +++ /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 25489) @@ -33,4 +33,42 @@ vector->Assemble(); + /*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; @@ -79,16 +117,33 @@ *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(); Index: /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h (revision 25488) +++ /issm/trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h (revision 25489) @@ -1,4 +1,4 @@ /*!\file: GetVectorFromInputsx.h - */ + */ #ifndef _GETVECTORFROMINPUTSXX_H @@ -8,8 +8,11 @@ /* 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 */