Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 22897) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp (revision 22898) @@ -100,4 +100,5 @@ iomodel->FetchDataToInput(elements,"md.initialization.sediment_head",SedimentHeadHydrostepEnum); iomodel->FetchDataToInput(elements,"md.hydrology.sediment_transmitivity",HydrologydcSedimentTransmitivityEnum); + iomodel->FetchDataToInput(elements,"md.hydrology.mask_thawed_node",HydrologydcMaskThawedNodeEnum); if(iomodel->domaintype!=Domain2DhorizontalEnum){ iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum); @@ -184,4 +185,5 @@ /*Intermediaries*/ + bool thawed_element; int domaintype; Element* basalelement; @@ -199,4 +201,17 @@ default: _error_("mesh "<GetInput(HydrologydcMaskThawedEltEnum); _assert_(thawed_element_input); + thawed_element_input->GetInputValue(&thawed_element); + + /*Check that all nodes are active, else return empty matrix*/ + if(!thawed_element) { + if(domaintype!=Domain2DhorizontalEnum){ + basalelement->DeleteMaterials(); + delete basalelement; + } + return NULL; + } + /*Intermediaries */ @@ -297,5 +312,6 @@ /*Intermediaries*/ - int domaintype; + bool thawed_element; + int domaintype; Element* basalelement; @@ -311,4 +327,16 @@ break; default: _error_("mesh "<GetInput(HydrologydcMaskThawedEltEnum); _assert_(thawed_element_input); + thawed_element_input->GetInputValue(&thawed_element); + + /*Check that all nodes are active, else return empty matrix*/ + if(!thawed_element) { + if(domaintype!=Domain2DhorizontalEnum){ + basalelement->DeleteMaterials(); + delete basalelement; + } + return NULL; } @@ -626,12 +654,4 @@ _error_("UnconfinedFlag is 0 or 1"); } - - /*Let's deal with the frozen parts for which transmitivity is zero*/ - if (isthermal){ - element->GetInputValue(&meltingrate,gauss,BasalforcingsGroundediceMeltingRateEnum); - element->GetInputValue(&groundedice,gauss,MaskGroundediceLevelsetEnum); - if ((meltingrate<=0.0) && (groundedice>0)) sediment_transmitivity=0.; - } - return sediment_transmitivity; }/*}}}*/ @@ -738,2 +758,63 @@ } }/*}}}*/ + +void HydrologyDCInefficientAnalysis::HydrologyIDSGetMask(Vector* vec_mask, Element* element){ + bool active_element; + int domaintype; + Element* basalelement=NULL; + + /*Get basal element*/ + element->FindParam(&domaintype,DomainTypeEnum); + switch(domaintype){ + case Domain2DhorizontalEnum: + basalelement = element; + break; + case Domain3DEnum: + if(!element->IsOnBase()) return; + basalelement = element->SpawnBasalElement(); + break; + default: _error_("mesh "<GetNumberOfNodes(); + IssmDouble* meltingrate = xNew(numnodes); + IssmDouble* groundedice = xNew(numnodes); + + basalelement->GetInputListOnVertices(&meltingrate[0],BasalforcingsGroundediceMeltingRateEnum); + basalelement->GetInputListOnVertices(&groundedice[0],MaskGroundediceLevelsetEnum); + + /*if melting rate is not positive and node is not floating, deactivate*/ + for(int i=0;i0)){ + vec_mask->SetValue(basalelement->nodes[i]->Sid(),0.,INS_VAL); + } + else{ + vec_mask->SetValue(basalelement->nodes[i]->Sid(),1.,INS_VAL); + } + } + + if(domaintype!=Domain2DhorizontalEnum){ + basalelement->DeleteMaterials(); + delete basalelement; + } + xDelete(meltingrate); + xDelete(groundedice); +} +void HydrologyDCInefficientAnalysis::ElementizeIdsMask(FemModel* femmodel){/*{{{*/ + + bool element_active; + Element* element=NULL; + + for(int i=0;ielements->Size();i++){ + element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + Input* node_mask_input = element->GetInput(HydrologydcMaskThawedNodeEnum); _assert_(node_mask_input); + + if(node_mask_input->Max()>0.){ + element_active = true; + } + else{ + element_active = false; + } + element->AddInput(new BoolInput(HydrologydcMaskThawedEltEnum,element_active)); + } +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h (revision 22897) +++ /issm/trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.h (revision 22898) @@ -1,3 +1,3 @@ -/*! \file HydrologyDCInefficientAnalysis.h +/*! \file HydrologyDCInefficientAnalysis.h * \brief: header file for generic external result object */ @@ -8,5 +8,5 @@ /*Headers*/ #include "./Analysis.h" -class Node; +class Node; class Input; class HydrologyDCInefficientAnalysis: public Analysis{ @@ -40,4 +40,6 @@ IssmDouble GetHydrologyPVectorTransfer(Element* element, Gauss* gauss, Input* epl_head_input); void ElementizeEplMask(FemModel* femmodel); + void HydrologyIDSGetMask(Vector* vec_mask, Element* element); + void ElementizeIdsMask(FemModel* femmodel); }; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22897) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 22898) @@ -2069,4 +2069,5 @@ name==HydrologydcEplThicknessHydrostepEnum || name==HydrologydcMaskEplactiveNodeEnum || + name==HydrologydcMaskThawedNodeEnum || name==HydrologyHeadEnum || name==HydrologyHeadOldEnum || Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22897) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 22898) @@ -136,6 +136,6 @@ bool traceon=true; this->profiler=NULL; /*avoid leak, as we are not using the profiler ever in ad control run. */ - - /*Store the communicator, but do not set it as a global variable, as this has already + + /*Store the communicator, but do not set it as a global variable, as this has already * been done by the FemModel that called this copy constructor: */ IssmComm::SetComm(incomm); @@ -152,5 +152,5 @@ this->amrbamg = NULL; #endif - + /*Save communicator in the parameters dataset: */ this->parameters->AddObject(new GenericParam(incomm,FemModelCommEnum)); @@ -182,5 +182,5 @@ if(parameters)delete parameters; if(results)delete results; - + #if defined(_HAVE_NEOPZ_) && !defined(_HAVE_ADOLC_) if(amr)delete amr; @@ -206,5 +206,5 @@ /*First, recover the name of the restart file: */ parameters->FindParam(&restartfilename,RestartFileNameEnum); - + /*Open file for writing: */ restartfid=pfopen(restartfilename,"wb"); @@ -215,8 +215,8 @@ /*Create buffer to hold marshalled femmodel: */ this->Marshall(NULL,&femmodel_size,MARSHALLING_SIZE); - femmodel_buffer=xNew(femmodel_size); + femmodel_buffer=xNew(femmodel_size); /*Keep track of initial position of femmodel_buffer: */ femmodel_buffer_ini=femmodel_buffer; - + /*Marshall:*/ this->Marshall(&femmodel_buffer,NULL,MARSHALLING_FORWARD); @@ -376,10 +376,10 @@ }/*}}}*/ void FemModel::InitFromFids(char* rootpath, FILE* IOMODEL, FILE* toolkitsoptionsfid, int in_solution_type, bool trace, IssmPDouble* X){/*{{{*/ - + /*Initialize internal data: */ this->solution_type = in_solution_type; this->analysis_counter = nummodels-1; //point to last analysis_type carried out. this->results = new Results(); //not initialized by CreateDataSets - + /*create IoModel */ IoModel* iomodel = new IoModel(IOMODEL,in_solution_type,trace,X); @@ -399,8 +399,8 @@ if(VerboseMProcessor()) _printf0_(" configuring element and loads\n"); ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters); - + if(i==0){ if(VerboseMProcessor()) _printf0_(" creating vertex PIDs\n"); - VerticesDofx(vertices,parameters); + VerticesDofx(vertices,parameters); if(VerboseMProcessor()) _printf0_(" detecting active vertices\n"); @@ -409,5 +409,5 @@ if(VerboseMProcessor()) _printf0_(" resolving node constraints\n"); - SpcNodesx(nodes,constraints,parameters,analysis_type_list[i]); + SpcNodesx(nodes,constraints,parameters,analysis_type_list[i]); if(VerboseMProcessor()) _printf0_(" creating nodal degrees of freedom\n"); @@ -487,6 +487,6 @@ FILE* restartfid=NULL; char* restartfilename = NULL; - int femmodel_size=0; - int fread_return=0; + int femmodel_size=0; + int fread_return=0; char* femmodel_buffer=NULL; char* femmodel_buffer_ini=NULL; @@ -513,10 +513,10 @@ /*Figure out size of buffer to be read: */ - fseek(restartfid, 0L, SEEK_END); + fseek(restartfid, 0L, SEEK_END); femmodel_size = ftell(restartfid); fseek(restartfid, 0L, SEEK_SET); /*Allocate buffer: */ - femmodel_buffer=xNew(femmodel_size); + femmodel_buffer=xNew(femmodel_size); /*Read buffer from file: */ @@ -539,7 +539,7 @@ void FemModel::SetCurrentConfiguration(int configuration_type,int analysis_type){/*{{{*/ - /*Use configuration_type to setup the analysis counter, the configurations of objects etc ... but use - * analysis_type to drive the element numerics. This allows for use of 1 configuration_type for several - * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the + /*Use configuration_type to setup the analysis counter, the configurations of objects etc ... but use + * analysis_type to drive the element numerics. This allows for use of 1 configuration_type for several + * analyses. For example: do a SurfaceSlopeX, SurfaceSlopeY, BedSlopeX and BedSlopeY analysis using the * Slope configuration.*/ int found=-1; @@ -700,5 +700,5 @@ analyses_temp[numanalyses++]=EsaAnalysisEnum; break; - + case SealevelriseSolutionEnum: analyses_temp[numanalyses++]=SealevelriseAnalysisEnum; @@ -825,11 +825,11 @@ /*run solution core: */ - profiler->Start(CORE); - solutioncore(this); + profiler->Start(CORE); + solutioncore(this); profiler->Stop(CORE); /*run AD core if needed: */ - profiler->Start(ADCORE); - ad_core(this); + profiler->Start(ADCORE); + ad_core(this); profiler->Stop(ADCORE); @@ -1121,5 +1121,5 @@ /*Broadcast and plug into response: */ ISSM_MPI_Allreduce ( &cpu_found,&cpu_found,1,ISSM_MPI_INT,ISSM_MPI_MAX,IssmComm::GetComm()); - ISSM_MPI_Bcast(&response,1,ISSM_MPI_DOUBLE,cpu_found,IssmComm::GetComm()); + ISSM_MPI_Bcast(&response,1,ISSM_MPI_DOUBLE,cpu_found,IssmComm::GetComm()); /*Assign output pointers: */ @@ -1268,5 +1268,5 @@ num_segments = mdims_array[counter-1]; - /*Go through segments, and then elements, and figure out which elements belong to a segment. + /*Go through segments, and then elements, and figure out which elements belong to a segment. * When we find one, use the element to compute the mass flux on the segment: */ for(i=0;iGetInputDerivativeValue(&dp[0],xyz_list,gauss); - /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ + /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ //J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight; J+=weight*1/2*pow(dp[0]*dp[0]+dp[1]*dp[1],2)*Jdet*gauss->weight; @@ -1678,5 +1678,5 @@ omega0_input->GetInputValue(&p0,gauss); - /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ + /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ //J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight; J+=weight*1/2*pow(p - p0,2)*Jdet*gauss->weight; @@ -1902,5 +1902,5 @@ ISSM_MPI_Bcast(&nodesperelement,1,ISSM_MPI_INT,0,IssmComm::GetComm()); ISSM_MPI_Bcast(&array_size,1,ISSM_MPI_INT,0,IssmComm::GetComm()); - + if(results_on_nodes){ @@ -1955,5 +1955,5 @@ IssmDouble* values = xNewZeroInit(nlines*ncols); IssmDouble* allvalues = xNew(nlines*ncols); - + /*Fill-in matrix*/ for(int j=0;jSize();j++){ @@ -1964,5 +1964,5 @@ ISSM_MPI_Allreduce((void*)values,(void*)allvalues,ncols*nlines,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); xDelete(values); - + if(save_results)results->AddResult(new GenericExternalResult(results->Size()+1,output_enum,allvalues,nlines,ncols,step,time)); xDelete(allvalues); @@ -2012,5 +2012,5 @@ int domaintype; this->parameters->FindParam(&domaintype,DomainTypeEnum); - + /*1: go throug all elements of this partition and figure out how many * segments we have (corresopnding to levelset = 0)*/ @@ -2132,11 +2132,11 @@ case VelEnum: this->ElementResponsex(responses,VelEnum); break; case FrictionCoefficientEnum: NodalValuex(responses, FrictionCoefficientEnum,elements,nodes, vertices, loads, materials, parameters); break; - default: + default: if(response_descriptor_enum>=Outputdefinition1Enum && response_descriptor_enum <=Outputdefinition100Enum){ IssmDouble double_result = OutputDefinitionsResponsex(this,response_descriptor_enum); *responses=double_result; } - else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); - break; + else _error_("response descriptor \"" << EnumToStringx(response_descriptor_enum) << "\" not supported yet!"); + break; } @@ -2269,5 +2269,5 @@ thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss); - /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ + /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ J+=weight*1/2*(dp[0]*dp[0]+dp[1]*dp[1])*Jdet*gauss->weight; } @@ -2460,7 +2460,7 @@ analysis->UpdateConstraints(this); delete analysis; - + /*Second, constraints might be time dependent: */ - SpcNodesx(nodes,constraints,parameters,analysis_type); + SpcNodesx(nodes,constraints,parameters,analysis_type); /*Now, update degrees of freedoms: */ @@ -2473,5 +2473,5 @@ IssmDouble *surface = NULL; IssmDouble *bed = NULL; - + if(VerboseSolution()) _printf0_(" updating vertices positions\n"); @@ -2512,5 +2512,5 @@ /*AMR*/ -#if !defined(_HAVE_ADOLC_) +#if !defined(_HAVE_ADOLC_) void FemModel::ReMesh(void){/*{{{*/ @@ -2522,5 +2522,5 @@ int newnumberofvertices = -1; int newnumberofelements = -1; - bool* my_elements = NULL; + bool* my_elements = NULL; int* my_vertices = NULL; int elementswidth = this->GetElementsWidth();//just tria elements in this version @@ -2528,5 +2528,5 @@ bool isgroundingline; - /*Branch to specific amr depending on requested method*/ + /*Branch to specific amr depending on requested method*/ this->parameters->FindParam(&amrtype,AmrTypeEnum); switch(amrtype){ @@ -2541,5 +2541,5 @@ default: _error_("not implemented yet"); } - + /*Partitioning the new mesh. Maybe ElementsAndVerticesPartitioning.cpp could be modified to set this without iomodel.*/ this->ElementsAndVerticesPartitioning(newnumberofvertices,newnumberofelements,elementswidth,newelementslist,&my_elements,&my_vertices); @@ -2572,5 +2572,5 @@ /*As the domain is 2D, it is not necessary to create nodes for this analysis*/ - if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue; + if(analysis_enum==StressbalanceVerticalAnalysisEnum) continue; this->CreateNodes(newnumberofvertices,my_vertices,nodecounter,analysis_enum,new_nodes); @@ -2602,5 +2602,5 @@ //SetCurrentConfiguration(analysis_type); - this->analysis_counter=i; + this->analysis_counter=i; /*Now, plug analysis_counter and analysis_type inside the parameters: */ this->parameters->SetParam(this->analysis_counter,AnalysisCounterEnum); @@ -2619,5 +2619,5 @@ ConfigureObjectsx(new_elements,this->loads,new_nodes,new_vertices,new_materials,this->parameters); - if(i==0){ + if(i==0){ VerticesDofx(new_vertices,this->parameters); //only call once, we only have one set of vertices } @@ -2663,5 +2663,5 @@ /*Insert bedrock from mismip+ setup*/ /*This was used to Misomip project/simulations*/ - + if(VerboseSolution())_printf0_(" call Mismip bedrock adjust module\n"); @@ -2680,5 +2680,5 @@ by = 500./(1.+exp((-2./4000.)*(y-80000./2.-24000.)))+500./(1.+exp((2./4000.)*(y-80000./2.+24000.))); r[i] = max(bx+by,-720.); - } + } /*insert new bedrock*/ element->AddInput(BedEnum,&r[0],P1Enum); @@ -2693,5 +2693,5 @@ if(VerboseSolution())_printf0_(" call adjust base and thickness module\n"); - + int numvertices = this->GetElementsWidth(); IssmDouble rho_water,rho_ice,density,base_float; @@ -2705,5 +2705,5 @@ for(int el=0;elelements->Size();el++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(el)); - + element->GetInputListOnVertices(&s[0],SurfaceEnum); element->GetInputListOnVertices(&r[0],BedEnum); @@ -2717,14 +2717,14 @@ base_float = rho_ice*s[i]/(rho_ice-rho_water); if(r[i]>base_float){ - b[i] = r[i]; - } + b[i] = r[i]; + } else { - b[i] = base_float; - } + b[i] = base_float; + } if(fabs(sl[i])>0) _error_("Sea level value "<AddInput(BaseEnum,&b[0],P1Enum); } - + /*Delete*/ xDelete(phi); @@ -2762,5 +2762,5 @@ Vector* input_interpolations = NULL; IssmDouble* input_interpolations_serial = NULL; - int* pos = NULL; + int* pos = NULL; IssmDouble value = 0; @@ -2782,5 +2782,5 @@ P0input_interp = xNew(numP0inputs); P1input_enums = xNew(numP1inputs); - P1input_interp = xNew(numP1inputs); + P1input_interp = xNew(numP1inputs); } numP0inputs = 0; @@ -2814,6 +2814,6 @@ } } - - /*Get P0 and P1 inputs over the elements*/ + + /*Get P0 and P1 inputs over the elements*/ pos = xNew(elementswidth); vP0inputs= new Vector(numberofelements*numP0inputs); @@ -2821,14 +2821,14 @@ for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); - + /*Get P0 inputs*/ for(int j=0;j(element->GetInput(P0input_enums[j])); + TriaInput* input=xDynamicCast(element->GetInput(P0input_enums[j])); input->GetInputAverage(&value); pos[0]=element->Sid()*numP0inputs+j; - /*Insert input in the vector*/ + /*Insert input in the vector*/ vP0inputs->SetValues(1,pos,&value,INS_VAL); } - + /*Get P1 inputs*/ for(int j=0;jvertices[1]->Sid()*numP1inputs+j; pos[2]=element->vertices[2]->Sid()*numP1inputs+j; - /*Insert input in the vector*/ - vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL); + /*Insert input in the vector*/ + vP1inputs->SetValues(elementswidth,pos,input->values,INS_VAL); } } @@ -2847,14 +2847,14 @@ P0inputs=vP0inputs->ToMPISerial(); P1inputs=vP1inputs->ToMPISerial(); - + /*Assign pointers*/ - *pnumP0inputs = numP0inputs; - *pP0inputs = P0inputs; + *pnumP0inputs = numP0inputs; + *pP0inputs = P0inputs; *pP0input_enums = P0input_enums; *pP0input_interp = P0input_interp; - *pnumP1inputs = numP1inputs; - *pP1inputs = P1inputs; + *pnumP1inputs = numP1inputs; + *pP1inputs = P1inputs; *pP1input_enums = P1input_enums; - *pP1input_interp = P1input_interp; + *pP1input_interp = P1input_interp; /*Cleanup*/ @@ -2867,5 +2867,5 @@ /*}}}*/ void FemModel::InterpolateInputs(Vertices* newfemmodel_vertices,Elements* newfemmodel_elements){/*{{{*/ - + int numberofelements = this->elements->NumberOfElements(); //global, entire old mesh int newnumberofelements = newfemmodel_elements->Size(); //just on the new partition @@ -2883,5 +2883,5 @@ int* P1input_enums = NULL; int* P1input_interp = NULL; - IssmDouble* values = NULL; + IssmDouble* values = NULL; IssmDouble* vector = NULL; IssmDouble* x = NULL;//global, entire old mesh @@ -2895,5 +2895,5 @@ IssmDouble* newyc = NULL;//just on the new partition int* newelementslist = NULL;//just on the new partition - int* sidtoindex = NULL;//global vertices sid to partition index + int* sidtoindex = NULL;//global vertices sid to partition index /*Get old P0 and P1 inputs (entire mesh)*/ @@ -2928,15 +2928,15 @@ /*Insert P0 and P1 inputs into the new elements (just on the new partition)*/ - values=xNew(elementswidth); + values=xNew(elementswidth); for(int i=0;iSize();i++){//just on the new partition Element* element=xDynamicCast(newfemmodel_elements->GetObjectByOffset(i)); /*newP0inputs is just on the new partition*/ for(int j=0;jAddInput(new DoubleInput(P0input_enums[j],newP0inputs[i*numP0inputs+j])); break; - case IntInputEnum: + case IntInputEnum: element->AddInput(new IntInput(P0input_enums[j],reCast(newP0inputs[i*numP0inputs+j]))); break; @@ -2956,5 +2956,5 @@ } } - + /*Cleanup*/ xDelete(P0inputs); @@ -2995,5 +2995,5 @@ if(!this->elements || !this->vertices || !this->results || !this->parameters) return; - + parameters->FindParam(&step,StepEnum); parameters->FindParam(&time,TimeEnum); @@ -3013,5 +3013,5 @@ this->results->AddResult(new GenericExternalResult(this->results->Size()+1,MeshYEnum, y,numberofvertices,1,step,time)); - + /*Cleanup*/ xDelete(x); @@ -3074,5 +3074,5 @@ /*Assemble*/ vmasklevelset->Assemble(); - + /*Serialize and set output*/ (*pmasklevelset)=vmasklevelset->ToMPISerial(); @@ -3088,5 +3088,5 @@ /*newelementslist is in Matlab indexing*/ - + /*Creating connectivity table*/ int* connectivity=NULL; @@ -3099,10 +3099,10 @@ connectivity[vertexid-1]+=1;//Matlab to C indexing } - } + } /*Create vertex and insert in vertices*/ for(int i=0;iid=i+1; newvertex->sid=i; @@ -3115,6 +3115,6 @@ newvertex->connectivity=connectivity[i]; newvertex->clone=false;//itapopo check this - vertices->AddObject(newvertex); - } + vertices->AddObject(newvertex); + } } @@ -3143,5 +3143,5 @@ } else newtria->element_type_list=NULL; - + /*Element hook*/ int matpar_id=newnumberofelements+1; //retrieve material parameter id (last pointer in femodel->materials) @@ -3149,5 +3149,5 @@ /*retrieve vertices ids*/ int* vertex_ids=xNew(elementswidth); - for(int j=0;j(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing + for(int j=0;j(newelementslist[elementswidth*i+j]);//this Hook wants Matlab indexing /*Setting the hooks*/ newtria->numanalyses =this->nummodels; @@ -3161,6 +3161,6 @@ /*Clean up*/ xDelete(vertex_ids); - elements->AddObject(newtria); - } + elements->AddObject(newtria); + } } @@ -3172,12 +3172,12 @@ for(int i=0;iAddObject(new Matice(i+1,i,MaticeEnum)); - } - } - + materials->AddObject(new Matice(i+1,i,MaticeEnum)); + } + } + /*Add new constant material property to materials, at the end: */ Matpar *newmatpar=static_cast(this->materials->GetObjectByOffset(this->materials->Size()-1)->copy()); newmatpar->SetMid(newnumberofelements+1); - materials->AddObject(newmatpar);//put it at the end of the materials + materials->AddObject(newmatpar);//put it at the end of the materials } /*}}}*/ @@ -3186,8 +3186,8 @@ int lid=0; for(int j=0;jid=nodecounter+j+1; @@ -3195,12 +3195,12 @@ newnode->lid=lid++; newnode->analysis_enum=analysis_enum; - + /*Initialize coord_system: Identity matrix by default*/ for(int k=0;k<3;k++) for(int l=0;l<3;l++) newnode->coord_system[k][l]=0.0; for(int k=0;k<3;k++) newnode->coord_system[k][k]=1.0; - + /*indexing:*/ newnode->indexingupdate=true; - + Analysis* analysis=EnumToAnalysis(analysis_enum); int *doftypes=NULL; @@ -3216,5 +3216,5 @@ /*Stressbalance Horiz*/ if(analysis_enum==StressbalanceAnalysisEnum){ - // itapopo this code is rarely used. + // itapopo this code is rarely used. /*Coordinate system provided, convert to coord_system matrix*/ //XZvectorsToCoordinateSystem(&this->coord_system[0][0],&iomodel->Data(StressbalanceReferentialEnum)[j*6]); @@ -3231,5 +3231,5 @@ if(!femmodel_vertices) _error_("GetMesh: vertices are NULL."); if(!femmodel_elements) _error_("GetMesh: elements are NULL."); - + int numberofvertices = femmodel_vertices->NumberOfVertices(); int numberofelements = femmodel_elements->NumberOfElements(); @@ -3237,5 +3237,5 @@ IssmDouble* x = NULL; IssmDouble* y = NULL; - IssmDouble* z = NULL; + IssmDouble* z = NULL; int* elementslist = NULL; int* elem_vertices = NULL; @@ -3246,5 +3246,5 @@ /*Get vertices coordinates*/ VertexCoordinatesx(&x,&y,&z,femmodel_vertices,false) ; - + /*Get element vertices*/ elem_vertices = xNew(elementswidth); @@ -3261,5 +3261,5 @@ vid3->SetValue(element->sid,elem_vertices[2],INS_VAL); } - + /*Assemble*/ vid1->Assemble(); @@ -3271,5 +3271,5 @@ id2 = vid2->ToMPISerial(); id3 = vid3->ToMPISerial(); - + /*Construct elements list*/ elementslist=xNew(numberofelements*elementswidth); @@ -3280,5 +3280,5 @@ elementslist[elementswidth*i+2] = reCast(id3[i])+1; //InterpMesh wants Matlab indexing } - + /*Assign pointers*/ *px = x; @@ -3301,5 +3301,5 @@ if(!femmodel_vertices) _error_("GetMesh: vertices are NULL."); if(!femmodel_elements) _error_("GetMesh: elements are NULL."); - + int numberofvertices = femmodel_vertices->Size(); //number of vertices of this partition int numbertotalofvertices = femmodel_vertices->NumberOfVertices(); //number total of vertices (entire mesh) @@ -3308,9 +3308,9 @@ IssmDouble* x = NULL; IssmDouble* y = NULL; - IssmDouble* z = NULL; + IssmDouble* z = NULL; int* elementslist = NULL; int* sidtoindex = NULL; int* elem_vertices = NULL; - + /*Get vertices coordinates of this partition*/ sidtoindex = xNewZeroInit(numbertotalofvertices);//entire mesh, all vertices @@ -3318,8 +3318,8 @@ y = xNew(numberofvertices);//just this partitio; z = xNew(numberofvertices);//just this partitio; - + /*Go through in this partition (vertices)*/ for(int i=0;iGetObjectByOffset(i); + Vertex* vertex=(Vertex*)femmodel_vertices->GetObjectByOffset(i); /*Attention: no spherical coordinates*/ x[i]=vertex->GetX(); @@ -3334,5 +3334,5 @@ elementslist = xNew(numberofelements*elementswidth); if(numberofelements*elementswidth<0) _error_("numberofelements negative."); - + for(int i=0;i(femmodel_elements->GetObjectByOffset(i)); @@ -3341,6 +3341,6 @@ elementslist[elementswidth*i+1] = sidtoindex[elem_vertices[1]]+1; //InterpMesh wants Matlab indexing elementslist[elementswidth*i+2] = sidtoindex[elem_vertices[2]]+1; //InterpMesh wants Matlab indexing - } - + } + /*Assign pointers*/ *px = x; @@ -3348,5 +3348,5 @@ *pz = z; *pelementslist = elementslist; //Matlab indexing. InterMesh uses this type. - *psidtoindex = sidtoindex; //it is ncessary to insert inputs + *psidtoindex = sidtoindex; //it is ncessary to insert inputs /*Cleanup*/ @@ -3359,7 +3359,7 @@ /*OTHERS CONSTRAINTS MUST BE IMPLEMENTED*/ if(analysis_enum!=StressbalanceAnalysisEnum) return; - + int numberofnodes_analysistype= this->nodes->NumberOfNodes(analysis_enum); - int dofpernode = 2; //vx and vy + int dofpernode = 2; //vx and vy int numberofcols = dofpernode*2; //to keep dofs and flags in the vspc vector int numberofvertices = this->vertices->NumberOfVertices(); //global, entire old mesh @@ -3385,5 +3385,5 @@ newy=xNew(newnumberofvertices);//just the new partition for(int i=0;iGetObjectByOffset(i); + Vertex* vertex=(Vertex*)newfemmodel_vertices->GetObjectByOffset(i); /*Attention: no spherical coordinates*/ newx[i]=vertex->GetX(); @@ -3393,5 +3393,5 @@ /*Get spcvx and spcvy of old mesh*/ for(int i=0;iconstraints->Size();i++){ - + Constraint* constraint=(Constraint*)constraints->GetObjectByOffset(i); if(!constraint->InAnalysis(analysis_enum)) _error_("AMR create constraints for "<GetDof(); int node = spcstatic->GetNodeId(); - IssmDouble spcvalue = spcstatic->GetValue(); + IssmDouble spcvalue = spcstatic->GetValue(); int nodeindex = node-1; - + /*vx and vx flag insertion*/ if(dof==0) {//vx vspc->SetValue(nodeindex*numberofcols,spcvalue,INS_VAL); //vx vspc->SetValue(nodeindex*numberofcols+dofpernode,1,INS_VAL);//vxflag - } + } /*vy and vy flag insertion*/ if(dof==1){//vy @@ -3423,5 +3423,5 @@ spc,numberofvertices,numberofcols, newx,newy,newnumberofvertices,NULL); - + /*Now, insert the interpolated constraints in the data set (constraints)*/ count=0; @@ -3440,5 +3440,5 @@ /*spcvy*/ if(!xIsNan(newspc[i*numberofcols+1]) && newspc[i*numberofcols+dofpernode+1]>(1-eps) ){ - newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum)); + newfemmodel_constraints->AddObject(new SpcStatic(constraintcounter+count+1,nodecounter+vertex->Sid()+1,1,newspc[i*numberofcols+1],analysis_enum)); //add count'th spc, on node i+1, setting dof 1 to vx. count++; @@ -3497,11 +3497,11 @@ bool *my_elements = NULL; int *my_vertices = NULL; - - _assert_(newnumberofvertices>0); - _assert_(newnumberofelements>0); + + _assert_(newnumberofvertices>0); + _assert_(newnumberofelements>0); epart=xNew(newnumberofelements); npart=xNew(newnumberofvertices); index=xNew(elementswidth*newnumberofelements); - + for (int i=0;i(newnumberofvertices); @@ -3533,12 +3533,12 @@ for(int i=0;i(epart); - xDelete(npart); + xDelete(npart); xDelete(index); } /*}}}*/ void FemModel::SmoothedDeviatoricStressTensor(IssmDouble** ptauxx,IssmDouble** ptauyy,IssmDouble** ptauxy){/*{{{*/ - + int elementswidth = this->GetElementsWidth();//just 2D mesh, tria elements int numberofvertices = this->vertices->NumberOfVertices(); IssmDouble weight = 0.; - IssmDouble* tauxx = NULL; - IssmDouble* tauyy = NULL; - IssmDouble* tauxy = NULL; + IssmDouble* tauxx = NULL; + IssmDouble* tauyy = NULL; + IssmDouble* tauxy = NULL; IssmDouble* totalweight = NULL; IssmDouble* deviatoricstressxx = xNew(elementswidth); @@ -3573,8 +3573,8 @@ Vector* vectauxy = new Vector(numberofvertices); Vector* vectotalweight = new Vector(numberofvertices); - + /*Update the Deviatoric Stress tensor over the elements*/ this->DeviatoricStressx(); - + /*Calculate the Smoothed Deviatoric Stress tensor*/ for(int i=0;ielements->Size();i++){ @@ -3621,5 +3621,5 @@ /*Divide for the total weight*/ for(int i=0;i0); + _assert_(totalweight[i]>0); tauxx[i] = tauxx[i]/totalweight[i]; tauyy[i] = tauyy[i]/totalweight[i]; @@ -3646,5 +3646,5 @@ void FemModel::ZZErrorEstimator(IssmDouble** pelementerror){/*{{{*/ - /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor. + /*Compute the Zienkiewicz and Zhu (ZZ) error estimator for the deviatoric stress tensor. * Ref.: Zienkiewicz and Zhu, A Simple Error Estimator and Adaptive Procedure for Practical Engineering Analysis, Int. J. Numer. Meth. Eng, 1987*/ @@ -3666,5 +3666,5 @@ /*Get smoothed deviatoric stress tensor*/ this->SmoothedDeviatoricStressTensor(&smoothedtauxx,&smoothedtauyy,&smoothedtauxy); - + /*Integrate the error over elements*/ for(int i=0;ielements->Size();i++){ @@ -3674,5 +3674,5 @@ element->GetInputListOnVertices(tauxy,DeviatoricStressxyEnum); element->GetVerticesSidList(elem_vertices); - + /*Integrate*/ element->GetVerticesCoordinates(&xyz_list); @@ -3689,7 +3689,7 @@ ftxy+=(tauxy[n]-smoothedtauxy[elem_vertices[n]])*basis[n]; } - error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2 - } - /*Set the error in the global vector*/ + error+=Jdet*gauss->weight*( pow(ftxx,2)+pow(ftyy,2)+pow(ftxy,2) ); //e^2 + } + /*Set the error in the global vector*/ sid=element->Sid(); error = sqrt(error);//sqrt(e^2) @@ -3705,5 +3705,5 @@ /*Serialize and set output*/ (*pelementerror)=velementerror->ToMPISerial(); - + /*Cleanup*/ xDelete(smoothedtauxx); @@ -3749,5 +3749,5 @@ Tria* triaelement = xDynamicCast(element); weight = triaelement->GetArea();//the tria area is a choice for the weight - + /*dH/dx*/ vecdHdx->SetValue(elem_vertices[0],weight*GradH[0],ADD_VAL); @@ -3817,5 +3817,5 @@ /*Get smoothed deviatoric stress tensor*/ this->SmoothedGradThickness(&smoothed_dHdx,&smoothed_dHdy); - + /*Integrate the error over elements*/ for(int i=0;ielements->Size();i++){ @@ -3903,5 +3903,5 @@ IssmDouble* x = NULL; IssmDouble* y = NULL; - IssmDouble* z = NULL; + IssmDouble* z = NULL; IssmDouble* xyz_list = NULL; IssmDouble x1,y1,x2,y2,x3,y3; @@ -3912,5 +3912,5 @@ //element->GetVerticesSidList(elem_vertices); int sid = element->Sid(); - element->GetVerticesCoordinates(&xyz_list); + element->GetVerticesCoordinates(&xyz_list); x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; @@ -3945,9 +3945,9 @@ _error_("level set type not implemented yet!"); } - + /*Outputs*/ IssmDouble* zerolevelset_points = NULL; int npoints = 0; - + /*Intermediaries*/ int elementswidth = this->GetElementsWidth(); @@ -3962,5 +3962,5 @@ int count,sid; IssmDouble xc,yc,x1,y1,x2,y2,x3,y3; - + /*Use the element center coordinate if level set is zero (grounding line or ice front), otherwise set NAN*/ for(int i=0;ielements->Size();i++){ @@ -3969,13 +3969,13 @@ element->GetVerticesSidList(elem_vertices); sid= element->Sid(); - element->GetVerticesCoordinates(&xyz_list); + element->GetVerticesCoordinates(&xyz_list); x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; x3 = xyz_list[3*2+0];y3 = xyz_list[3*2+1]; xc = NAN; - yc = NAN; + yc = NAN; Tria* tria = xDynamicCast(element); if(tria->IsIceInElement()){/*verify if there is ice in the element*/ - if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. || + if(levelset[0]*levelset[1]<0. || levelset[0]*levelset[2]<0. || abs(levelset[0]*levelset[1])Size(); - + /*Figure out max of ns: */ ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm()); @@ -4162,5 +4162,5 @@ } } - + /*One last time: */ pUp->Assemble(); @@ -4181,9 +4181,9 @@ int ns,nsmax; - + /*Go through elements, and add contribution from each element to the deflection vector wg:*/ ns = elements->Size(); - - /*First, figure out the surface area of Earth: */ + + /*First, figure out the surface area of Earth: */ for(int i=0;i(elements->GetObjectByOffset(i)); @@ -4209,5 +4209,5 @@ } } - + /*One last time: */ pUp->Assemble(); @@ -4237,5 +4237,5 @@ IssmDouble eartharea_cpu = 0.; int ns,nsmax; - + /*Go through elements, and add contribution from each element to the deflection vector wg:*/ ns = elements->Size(); @@ -4261,7 +4261,7 @@ for(int i=0;i(elements->GetObjectByOffset(i)); element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea); @@ -4271,5 +4271,5 @@ } if(VerboseConvergence())_printf0_("\n"); - + /*One last time: */ pSgi->Assemble(); @@ -4289,10 +4289,10 @@ /*serialized vectors:*/ IssmDouble* Sg_old=NULL; - + IssmDouble eartharea=0; IssmDouble eartharea_cpu=0; int ns,nsmax; - + /*Serialize vectors from previous iteration:*/ Sg_old=pSg_old->ToMPISerial(); @@ -4306,5 +4306,5 @@ eartharea_cpu += element->GetAreaSpherical(); } - + ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); @@ -4324,5 +4324,5 @@ } if(verboseconvolution)if(VerboseConvergence())_printf_("\n"); - + /*Free ressources:*/ xDelete(Sg_old); @@ -4336,7 +4336,7 @@ IssmDouble eartharea_cpu=0; IssmDouble tide_love_h, tide_love_k, fluid_love, moi_e, moi_p, omega, g; - IssmDouble load_love_k2 = -0.30922675; //degree 2 load Love number - IssmDouble m1, m2, m3; - IssmDouble lati, longi, radi, value; + IssmDouble load_love_k2 = -0.30922675; //degree 2 load Love number + IssmDouble m1, m2, m3; + IssmDouble lati, longi, radi, value; /*Serialize vectors from previous iteration:*/ @@ -4351,22 +4351,22 @@ ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - IssmDouble moi_list[3]={0,0,0}; - IssmDouble moi_list_cpu[3]={0,0,0}; + IssmDouble moi_list[3]={0,0,0}; + IssmDouble moi_list_cpu[3]={0,0,0}; for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); element->SealevelriseMomentOfInertia(&moi_list[0],Sg_old,eartharea); - moi_list_cpu[0] += moi_list[0]; - moi_list_cpu[1] += moi_list[1]; - moi_list_cpu[2] += moi_list[2]; + moi_list_cpu[0] += moi_list[0]; + moi_list_cpu[1] += moi_list[1]; + moi_list_cpu[2] += moi_list[2]; } ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - // + // ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - // + // ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - + /*pull out some useful parameters: */ parameters->FindParam(&tide_love_h,SealevelriseTidalLoveHEnum); @@ -4378,12 +4378,12 @@ /*compute perturbation terms for angular velocity vector: */ - m1 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[0]; - m2 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[1]; - m3 = -(1+load_love_k2)/moi_p * moi_list[2]; // term associated with fluid number (3-order-of-magnitude smaller) is negelected + m1 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[0]; + m2 = 1/(1-tide_love_k/fluid_love) * (1+load_love_k2)/(moi_p-moi_e) * moi_list[1]; + m3 = -(1+load_love_k2)/moi_p * moi_list[2]; // term associated with fluid number (3-order-of-magnitude smaller) is negelected /* Green's function (1+k_2-h_2/g): checked against Glenn Milne's thesis Chapter 3 (eqs: 3.3-4, 3.10-11) - * Perturbation terms for angular velocity vector (m1, m2, m3): checked against Mitrovica (2005 Appendix) & Jensen et al (2013 Appendix A3) - * Sea level rotational feedback: checked against GMD eqs 8-9 (only first order terms, i.e., degree 2 order 0 & 1 considered) - * all DONE in Geographic coordinates: theta \in [-90,90], lambda \in [-180 180] + * Perturbation terms for angular velocity vector (m1, m2, m3): checked against Mitrovica (2005 Appendix) & Jensen et al (2013 Appendix A3) + * Sea level rotational feedback: checked against GMD eqs 8-9 (only first order terms, i.e., degree 2 order 0 & 1 considered) + * all DONE in Geographic coordinates: theta \in [-90,90], lambda \in [-180 180] */ for(int i=0;iSize();i++){ @@ -4395,8 +4395,8 @@ lati=latitude[sid]/180*PI; longi=longitude[sid]/180*PI; radi=radius[sid]; - /*only first order terms are considered now: */ + /*only first order terms are considered now: */ value=((1.0+tide_love_k-tide_love_h)/9.81)*pow(omega*radi,2.0)* - (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi))); - + (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi))); + pSgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times } @@ -4404,5 +4404,5 @@ /*Assemble mesh velocity*/ pSgo_rot->Assemble(); - + /*Assign output pointers:*/ *pIxz=moi_list[0]; @@ -4412,5 +4412,5 @@ /*Free ressources:*/ xDelete(Sg_old); - + } /*}}}*/ @@ -4419,10 +4419,10 @@ /*serialized vectors:*/ IssmDouble* Sg=NULL; - + IssmDouble eartharea=0; IssmDouble eartharea_cpu=0; int ns,nsmax; - + /*Serialize vectors from previous iteration:*/ Sg=pSg->ToMPISerial(); @@ -4430,5 +4430,5 @@ /*Go through elements, and add contribution from each element to the deflection vector wg:*/ ns = elements->Size(); - + /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */ for(int i=0;iAssemble(); @@ -4492,5 +4492,5 @@ ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - + ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); @@ -4498,5 +4498,5 @@ /*Free ressources:*/ xDelete(Sg_serial); - + return oceanvalue/oceanarea; } @@ -4514,5 +4514,5 @@ int* eplzigzag_counter = NULL; int eplflip_lock; - + HydrologyDCEfficientAnalysis* effanalysis = new HydrologyDCEfficientAnalysis(); HydrologyDCInefficientAnalysis* inefanalysis = new HydrologyDCInefficientAnalysis(); @@ -4521,8 +4521,8 @@ mask=new Vector(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); recurence=new Vector(this->nodes->NumberOfNodes(HydrologyDCEfficientAnalysisEnum)); - this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); - this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); + this->parameters->FindParam(&eplzigzag_counter,NULL,EplZigZagCounterEnum); + this->parameters->FindParam(&eplflip_lock,HydrologydcEplflipLockEnum); GetVectorFromInputsx(&old_active,this,HydrologydcMaskEplactiveNodeEnum,NodeSIdEnum); - + for (int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); @@ -4542,6 +4542,6 @@ /*Assemble and serialize*/ mask->Assemble(); - serial_mask=mask->ToMPISerial(); - + serial_mask=mask->ToMPISerial(); + xDelete(eplzigzag_counter); xDelete(serial_rec); @@ -4585,8 +4585,67 @@ int sum_counter; ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); counter=sum_counter; *pEplcount = counter; if(VerboseSolution()) _printf0_(" Number of active nodes in EPL layer: "<< counter <<"\n"); + + /*Update dof indexings*/ + this->UpdateConstraintsx(); + +} +/*}}}*/ +void FemModel::HydrologyIDSupdateDomainx(IssmDouble* pIDScount){ /*{{{*/ + + bool isthermal; + Vector* mask = NULL; + IssmDouble* serial_mask = NULL; + + HydrologyDCInefficientAnalysis* inefanalysis = new HydrologyDCInefficientAnalysis(); + parameters->FindParam(&isthermal,TransientIsthermalEnum); + + /*When solving a thermal model we update the thawed nodes*/ + if(isthermal){ + /*Step 1: update mask, the mask correspond to thawed nodes (that have a meltingrate)*/ + mask=new Vector(this->nodes->NumberOfNodes(HydrologyDCInefficientAnalysisEnum)); + + for (int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + inefanalysis->HydrologyIDSGetMask(mask,element); + } + /*Assemble and serialize*/ + mask->Assemble(); + serial_mask=mask->ToMPISerial(); + delete mask; + } + /*for other cases we just grab the mask from the initialisation value*/ + else{ + GetVectorFromInputsx(&serial_mask,this,HydrologydcMaskThawedNodeEnum,NodeSIdEnum); + } + /*Update node activation accordingly*/ + int counter =0; + for (int i=0;iSize();i++){ + Node* node=xDynamicCast(nodes->GetObjectByOffset(i)); + if(node->InAnalysis(HydrologyDCInefficientAnalysisEnum)){ + if(serial_mask[node->Sid()]==1.){ + node->Activate(); + if(!node->IsClone()) counter++; + } + else{ + node->Deactivate(); + } + } + } + + /*Update Mask*/ + InputUpdateFromVectorx(this,serial_mask,HydrologydcMaskThawedNodeEnum,NodeSIdEnum); + xDelete(serial_mask); + inefanalysis->ElementizeIdsMask(this); + + int sum_counter; + ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + counter=sum_counter; + *pIDScount = counter; + if(VerboseSolution()) _printf0_(" Number of active nodes in IDS layer: "<< counter <<"\n"); /*Update dof indexings*/ @@ -4631,5 +4690,5 @@ int sum_counter; ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(&sum_counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); counter=sum_counter; *pL2count = counter; @@ -4716,5 +4775,5 @@ } /*}}}*/ -#ifdef _HAVE_JAVASCRIPT_ +#ifdef _HAVE_JAVASCRIPT_ FemModel::FemModel(IssmDouble* buffer, int buffersize, char* toolkits, char* solution, char* modelname,ISSM_MPI_Comm incomm, bool trace){ /*{{{*/ /*configuration: */ @@ -4731,10 +4790,10 @@ /*From command line arguments, retrieve different filenames needed to create the FemModel: */ solution_type=StringToEnumx(solution); - + /*Create femmodel from input files: */ profiler->Start(MPROCESSOR); this->InitFromBuffers((char*)buffer,buffersize,toolkits, solution_type,trace,NULL); profiler->Stop(MPROCESSOR); - + /*Save communicator in the parameters dataset: */ this->parameters->AddObject(new GenericParam(incomm,FemModelCommEnum)); @@ -4751,5 +4810,5 @@ size_t* poutputbuffersize; - + /*Before we delete the profiler, report statistics for this run: */ profiler->Stop(TOTAL); //final tagging @@ -4764,5 +4823,5 @@ ); _printf0_("\n"); - + /*Before we close the output file, recover the buffer and size:*/ outputbufferparam = xDynamicCast*>(this->parameters->FindParamObject(OutputBufferPointerEnum)); @@ -4804,5 +4863,5 @@ /*Open output file once for all and add output file descriptor to parameters*/ - output_fid=open_memstream(&outputbuffer,&outputsize); + output_fid=open_memstream(&outputbuffer,&outputsize); if(output_fid==NULL)_error_("could not initialize output stream"); this->parameters->SetParam(output_fid,OutputFilePointerEnum); @@ -4845,5 +4904,5 @@ /*Initialize hmaxvertices with NAN*/ hmaxvertices_serial=xNew(numberofvertices); - for(int i=0;iamrbamg->groundingline_distance>0) this->GethmaxVerticesFromZeroLevelSetDistance(hmaxvertices_serial,MaskGroundediceLevelsetEnum); @@ -4852,5 +4911,5 @@ if(this->amrbamg->deviatoricerror_threshold>0) this->GethmaxVerticesFromEstimators(hmaxvertices_serial,DeviatoricStressErrorEstimatorEnum); } - + if(my_rank==0){ this->amrbamg->ExecuteRefinementBamg(vector_serial,hmaxvertices_serial,&newnumberofvertices,&newnumberofelements,&newx,&newy,&newz,&newelementslist); @@ -4872,8 +4931,8 @@ newelementslist=xNew(newnumberofelements*this->GetElementsWidth()); } - ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); /*Assign output pointers*/ @@ -4908,5 +4967,5 @@ /*Create bamg data structures for bamg*/ this->amrbamg = new AmrBamg(); - + /*Get amr parameters*/ this->parameters->FindParam(&hmin,AmrHminEnum); @@ -4932,5 +4991,5 @@ /*Re-create original mesh and put it in bamg structure (only cpu 0)*/ - if(my_rank==0){ + if(my_rank==0){ this->amrbamg->Initialize(elements,x,y,numberofvertices,numberofelements); } @@ -4946,5 +5005,5 @@ if(!hmaxvertices) _error_("hmaxvertices is NULL!\n"); - + /*Intermediaries*/ int numberofvertices = this->vertices->NumberOfVertices(); @@ -4953,5 +5012,5 @@ switch(levelset_type){ - case MaskGroundediceLevelsetEnum: + case MaskGroundediceLevelsetEnum: threshold = this->amrbamg->groundingline_distance; resolution = this->amrbamg->groundingline_resolution; @@ -4967,5 +5026,5 @@ this->GetVerticeDistanceToZeroLevelSet(&verticedistance,levelset_type); if(!verticedistance) _error_("verticedistance is NULL!\n"); - + /*Fill hmaxVertices*/ for(int i=0;ivertices->NumberOfVertices(); IssmDouble* maxlength = xNew(numberofelements); - IssmDouble* error_vertices = xNewZeroInit(numberofvertices); + IssmDouble* error_vertices = xNewZeroInit(numberofvertices); IssmDouble* error_elements = NULL; IssmDouble* x = NULL; @@ -5002,5 +5061,5 @@ /*Fill variables*/ switch(errorestimator_type){ - case ThicknessErrorEstimatorEnum: + case ThicknessErrorEstimatorEnum: threshold = this->amrbamg->thicknesserror_threshold; groupthreshold = this->amrbamg->thicknesserror_groupthreshold; @@ -5027,10 +5086,10 @@ case ThicknessErrorEstimatorEnum: this->amrbamg->thicknesserror_maximum = maxerror;break; case DeviatoricStressErrorEstimatorEnum: this->amrbamg->deviatoricerror_maximum = maxerror;break; - } + } } /*Get mesh*/ this->GetMesh(this->vertices,this->elements,&x,&y,&z,&index); - + /*Fill error_vertices (this is the sum of all elements connected to the vertex)*/ for(int i=0;ivertices->NumberOfVertices(); @@ -5102,6 +5161,6 @@ /*Get vertices coordinates*/ VertexCoordinatesx(&x,&y,&z,this->vertices,false) ; - - /*Get points which level set is zero (center of elements with zero level set)*/ + + /*Get points which level set is zero (center of elements with zero level set)*/ this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type); @@ -5112,7 +5171,7 @@ for(int j=0;jamr->groundingline_distance>0) this->GetElementDistanceToZeroLevelSet(&gl_distance,MaskGroundediceLevelsetEnum); if(this->amr->icefront_distance>0) this->GetElementDistanceToZeroLevelSet(&if_distance,MaskIceLevelsetEnum); - if(this->amr->thicknesserror_threshold>0) this->ThicknessZZErrorEstimator(&thicknesserror); - if(this->amr->deviatoricerror_threshold>0) this->ZZErrorEstimator(&deviatoricerror); - + if(this->amr->thicknesserror_threshold>0) this->ThicknessZZErrorEstimator(&thicknesserror); + if(this->amr->deviatoricerror_threshold>0) this->ZZErrorEstimator(&deviatoricerror); + if(my_rank==0){ this->amr->ExecuteRefinement(gl_distance,if_distance,deviatoricerror,thicknesserror, - &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); + &newnumberofvertices,&newnumberofelements,&newx,&newy,&newelementslist); newz=xNewZeroInit(newnumberofvertices); if(newnumberofvertices<=0 || newnumberofelements<=0) _error_("Error in the ReMeshNeopz."); @@ -5167,10 +5226,10 @@ newelementslist=xNew(newnumberofelements*this->GetElementsWidth()); } - ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); - - /*Assign the pointers*/ + ISSM_MPI_Bcast(newx,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newy,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newz,newnumberofvertices,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + ISSM_MPI_Bcast(newelementslist,newnumberofelements*this->GetElementsWidth(),ISSM_MPI_INT,0,IssmComm::GetComm()); + + /*Assign the pointers*/ (*pnewelementslist) = newelementslist; //Matlab indexing (*pnewx) = newx; @@ -5188,5 +5247,5 @@ /*}}}*/ void FemModel::InitializeAdaptiveRefinementNeopz(void){/*{{{*/ - + /*Define variables*/ int my_rank = IssmComm::GetRank(); @@ -5197,5 +5256,5 @@ IssmDouble* z = NULL; int* elements = NULL; - + /*Initialize field as NULL for now*/ this->amr = NULL; @@ -5205,10 +5264,10 @@ this->GetMesh(this->vertices,this->elements,&x,&y,&z,&elements); - /*Create initial mesh (coarse mesh) in neopz data structure*/ + /*Create initial mesh (coarse mesh) in neopz data structure*/ /*Just CPU #0 should keep AMR object*/ /*Initialize refinement pattern*/ this->SetRefPatterns(); this->amr = new AdaptiveMeshRefinement(); - this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster) + this->amr->refinement_type=1;//1 is refpattern; 0 is uniform (faster) /*Get amr parameters*/ this->parameters->FindParam(&this->amr->level_max,AmrLevelMaxEnum); @@ -5223,5 +5282,5 @@ this->parameters->FindParam(&this->amr->deviatoricerror_groupthreshold,AmrDeviatoricErrorGroupThresholdEnum); this->parameters->FindParam(&this->amr->deviatoricerror_maximum,AmrDeviatoricErrorMaximumEnum); - if(my_rank==0){ + if(my_rank==0){ this->amr->CreateInitialMesh(numberofvertices,numberofelements,x,y,elements); } @@ -5244,5 +5303,5 @@ /*Output*/ IssmDouble* elementdistance; - + /*Intermediaries*/ int numberofelements = this->elements->NumberOfElements(); @@ -5253,6 +5312,6 @@ IssmDouble xc,yc,x1,y1,x2,y2,x3,y3; int numberofpoints; - - /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/ + + /*Get points which level set is zero (center of elements with zero level set, levelset_points is serial)*/ this->GetZeroLevelSetPoints(&levelset_points,numberofpoints,levelset_type); @@ -5260,5 +5319,5 @@ Element* element=xDynamicCast(this->elements->GetObjectByOffset(i)); int sid = element->Sid(); - element->GetVerticesCoordinates(&xyz_list); + element->GetVerticesCoordinates(&xyz_list); x1 = xyz_list[3*0+0];y1 = xyz_list[3*0+1]; x2 = xyz_list[3*1+0];y2 = xyz_list[3*1+1]; @@ -5270,9 +5329,9 @@ for(int j=0;jSetValue(sid,mindistance,INS_VAL); xDelete(xyz_list); - } + } /*Assemble*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 22897) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 22898) @@ -1,4 +1,4 @@ /* - * FemModel.h: + * FemModel.h: */ @@ -47,5 +47,5 @@ Nodes *nodes; //one set of nodes Parameters *parameters; //one set of parameters, independent of the analysis_type - Results *results; //results that cannot be fit into the elements + Results *results; //results that cannot be fit into the elements Vertices *vertices; //one set of vertices @@ -80,5 +80,5 @@ void Solve(void); - /*Modules*/ + /*Modules*/ void BalancethicknessMisfitx(IssmDouble* pV); void CalvingRateVonmisesx(); @@ -135,6 +135,6 @@ #endif #ifdef _HAVE_ESA_ - void EsaGeodetic2D(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pX, Vector* pY, IssmDouble* xx, IssmDouble* yy); - void EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); + void EsaGeodetic2D(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pX, Vector* pY, IssmDouble* xx, IssmDouble* yy); + void EsaGeodetic3D(Vector* pUp, Vector* pNorth, Vector* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); #endif #ifdef _HAVE_SEALEVELRISE_ @@ -142,8 +142,9 @@ void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution); void SealevelriseRotationalFeedback(Vector* pSgo_rot, Vector* pSg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); - void SealevelriseGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); + void SealevelriseGeodetic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); IssmDouble SealevelriseOceanAverage(Vector* Sg); #endif void HydrologyEPLupdateDomainx(IssmDouble* pEplcount); + void HydrologyIDSupdateDomainx(IssmDouble* pIDScount); void TimeAdaptx(IssmDouble* pdt); void UpdateConstraintsExtrudeFromBasex(); Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 22897) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 22898) @@ -13,10 +13,11 @@ /*intermediary*/ - int hydrology_model; - int solution_type; - int numoutputs=0; - bool save_results; - bool modify_loads=true; - char **requested_outputs = NULL; + int hydrology_model; + int solution_type; + int numoutputs = 0; + bool save_results; + bool modify_loads = true; + char **requested_outputs = NULL; + IssmDouble ThawedNodes; /*first recover parameters common to all solutions*/ @@ -45,8 +46,9 @@ else if (hydrology_model==HydrologydcEnum){ /*intermediary: */ - bool isefficientlayer; - int step,hydroslices; - IssmDouble time,init_time,hydrotime,yts; - IssmDouble dt,hydrodt; + bool isefficientlayer; + bool isthermal; + int step,hydroslices; + IssmDouble time,init_time,hydrotime,yts; + IssmDouble dt,hydrodt; femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); @@ -57,58 +59,64 @@ femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); - init_time = time-dt; //getting the time back to the start of the timestep - hydrotime=init_time; - hydrodt=dt/hydroslices; //computing hydro dt from dt and a divider - femmodel->parameters->AddObject(new DoubleParam(HydrologydtEnum,hydrodt)); - if(hydroslices>1){ - /*define which variable needs to be averaged on the sub-timestep and initialize as needed*/ - if (isefficientlayer){ - int inputtostack[4]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum,EplHeadHydrostepEnum,HydrologydcEplThicknessHydrostepEnum}; - int stackedinput[4]={EffectivePressureStackedEnum,SedimentHeadStackedEnum,EplHeadStackedEnum,HydrologydcEplThicknessStackedEnum}; - int averagedinput[4]={EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum}; - femmodel->InitTransientOutputx(&stackedinput[0],4); - while(hydrotimeparameters->FindParam(&isthermal,TransientIsthermalEnum); + /*first we exclude frozen nodes of the solved nodes*/ + femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum); + femmodel->HydrologyIDSupdateDomainx(&ThawedNodes); + if(ThawedNodes>0){ + init_time = time-dt; //getting the time back to the start of the timestep + hydrotime=init_time; + hydrodt=dt/hydroslices; //computing hydro dt from dt and a divider + femmodel->parameters->AddObject(new DoubleParam(HydrologydtEnum,hydrodt)); + if(hydroslices>1){ + /*define which variable needs to be averaged on the sub-timestep and initialize as needed*/ + if (isefficientlayer){ + int inputtostack[4]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum,EplHeadHydrostepEnum,HydrologydcEplThicknessHydrostepEnum}; + int stackedinput[4]={EffectivePressureStackedEnum,SedimentHeadStackedEnum,EplHeadStackedEnum,HydrologydcEplThicknessStackedEnum}; + int averagedinput[4]={EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum}; + femmodel->InitTransientOutputx(&stackedinput[0],4); + while(hydrotimeStackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,4); + } + femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,4); + } + else{ + int inputtostack[2]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum}; + int stackedinput[2]={EffectivePressureStackedEnum,SedimentHeadStackedEnum}; + int averagedinput[2]={EffectivePressureEnum,SedimentHeadEnum}; + femmodel->InitTransientOutputx(&stackedinput[0],2); + while(hydrotimeStackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,2); + } + femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,2); + } + } + else{ + InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum); + if (isefficientlayer){ InputDuplicatex(femmodel,EplHeadHydrostepEnum,EplHeadOldEnum); InputDuplicatex(femmodel,HydrologydcEplThicknessHydrostepEnum,HydrologydcEplThicknessOldEnum); - /*Proceed now to heads computations*/ - solutionsequence_hydro_nonlinear(femmodel); - /*If we have a sub-timestep we stack the variables here*/ - femmodel->StackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,4); } - femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,4); + /*Proceed now to heads computations*/ + solutionsequence_hydro_nonlinear(femmodel); } - else{ - int inputtostack[2]={EffectivePressureHydrostepEnum,SedimentHeadHydrostepEnum}; - int stackedinput[2]={EffectivePressureStackedEnum,SedimentHeadStackedEnum}; - int averagedinput[2]={EffectivePressureEnum,SedimentHeadEnum}; - femmodel->InitTransientOutputx(&stackedinput[0],2); - while(hydrotimeStackTransientOutputx(&inputtostack[0],&stackedinput[0],hydrotime,2); - } - femmodel->AverageTransientOutputx(&stackedinput[0],&averagedinput[0],init_time,2); - } - } - else{ - InputDuplicatex(femmodel,SedimentHeadHydrostepEnum,SedimentHeadOldEnum); - if (isefficientlayer){ - InputDuplicatex(femmodel,EplHeadHydrostepEnum,EplHeadOldEnum); - InputDuplicatex(femmodel,HydrologydcEplThicknessHydrostepEnum,HydrologydcEplThicknessOldEnum); - } - /*Proceed now to heads computations*/ - solutionsequence_hydro_nonlinear(femmodel); } } else if (hydrology_model==HydrologysommersEnum){ femmodel->SetCurrentConfiguration(HydrologySommersAnalysisEnum); - InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum); + InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum); solutionsequence_shakti_nonlinear(femmodel); if(VerboseSolution()) _printf0_(" updating gap height\n"); @@ -122,5 +130,9 @@ if(save_results){ if(VerboseSolution()) _printf0_(" saving results \n"); - femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); + if(hydrology_model==HydrologydcEnum && ThawedNodes==0){ + if(VerboseSolution()) _printf0_(" No thawed node hydro is skiped \n");} + else{ + femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); + } } /*Free ressources:*/ Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 22897) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 22898) @@ -1,5 +1,5 @@ /*!\file: transient_3d_core.cpp - * \brief: core of the transient_3d solution - */ + * \brief: core of the transient_3d solution + */ #ifdef HAVE_CONFIG_H @@ -142,5 +142,5 @@ InterpFromMeshToMesh2dx(&base_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice, icebase,ngrids_ice,1, - oceangridx,oceangridy,ngrids_ocean,options); + oceangridx,oceangridy,ngrids_ocean,options); delete options; @@ -239,5 +239,5 @@ int ngrids_ice=femmodel->vertices->NumberOfVertices(); int nels_ice=femmodel->elements->NumberOfElements(); - + /*Recover mesh and inputs needed*/ femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum); @@ -253,6 +253,6 @@ InterpFromMeshToMesh2dx(&base_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice, icebase,ngrids_ice,1, - oceangridx,oceangridy,ngrids_ocean,NULL); - + oceangridx,oceangridy,ngrids_ocean,NULL); + /*Send and receive data*/ ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm); @@ -267,5 +267,5 @@ InterpFromMeshToMesh2dx(&melt_mesh,index_ocean,oceangridx,oceangridy,ngrids_ocean,nels_ocean, oceanmelt,ngrids_ocean,1, - lon_ice,lat_ice,ngrids_ice,NULL); + lon_ice,lat_ice,ngrids_ice,NULL); for(int i=0;iUpdateVertexPositionsx(); } - + if(isgroundingline){ if(VerboseSolution()) _printf0_(" computing new grounding line position\n"); @@ -338,5 +338,5 @@ femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum); extrudefrombase_core(femmodel); - + if(save_results){ int outputs[3] = {SurfaceEnum,BaseEnum,MaskGroundediceLevelsetEnum}; @@ -356,5 +356,5 @@ /*esa: */ if(isesa) esa_core(femmodel); - + /*Sea level rise: */ if(isslr || iscoupler) sealevelrise_core(femmodel); @@ -367,5 +367,5 @@ femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1,save_results); } - + if(save_results){ if(VerboseSolution()) _printf0_(" saving temporary results\n"); @@ -402,9 +402,9 @@ } } - + if(!iscontrol || !isautodiff) femmodel->RequestedDependentsx(); if(iscontrol && isautodiff) femmodel->parameters->SetParam(dependent_objects,AutodiffDependentObjectsEnum); - /*Free ressources:*/ + /*Free ressources:*/ if(numoutputs){for(int i=0;i(requested_outputs[i]);} xDelete(requested_outputs);} } Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22897) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 22898) @@ -441,4 +441,6 @@ HydrologydcMaskEplactiveEltEnum, HydrologydcMaskEplactiveNodeEnum, + HydrologydcMaskThawedEltEnum, + HydrologydcMaskThawedNodeEnum, HydrologydcSedimentTransmitivityEnum, HydrologyEnglacialInputEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22897) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 22898) @@ -447,4 +447,6 @@ case HydrologydcMaskEplactiveEltEnum : return "HydrologydcMaskEplactiveElt"; case HydrologydcMaskEplactiveNodeEnum : return "HydrologydcMaskEplactiveNode"; + case HydrologydcMaskThawedEltEnum : return "HydrologydcMaskThawedElt"; + case HydrologydcMaskThawedNodeEnum : return "HydrologydcMaskThawedNode"; case HydrologydcSedimentTransmitivityEnum : return "HydrologydcSedimentTransmitivity"; case HydrologyEnglacialInputEnum : return "HydrologyEnglacialInput"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22897) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 22898) @@ -456,4 +456,6 @@ else if (strcmp(name,"HydrologydcMaskEplactiveElt")==0) return HydrologydcMaskEplactiveEltEnum; else if (strcmp(name,"HydrologydcMaskEplactiveNode")==0) return HydrologydcMaskEplactiveNodeEnum; + else if (strcmp(name,"HydrologydcMaskThawedElt")==0) return HydrologydcMaskThawedEltEnum; + else if (strcmp(name,"HydrologydcMaskThawedNode")==0) return HydrologydcMaskThawedNodeEnum; else if (strcmp(name,"HydrologydcSedimentTransmitivity")==0) return HydrologydcSedimentTransmitivityEnum; else if (strcmp(name,"HydrologyEnglacialInput")==0) return HydrologyEnglacialInputEnum; @@ -504,10 +506,10 @@ else if (strcmp(name,"Pressure")==0) return PressureEnum; else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum; - else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; - else if (strcmp(name,"Sealevel")==0) return SealevelEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum; + if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; + else if (strcmp(name,"Sealevel")==0) return SealevelEnum; + else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum; else if (strcmp(name,"SedimentHeadHydrostep")==0) return SedimentHeadHydrostepEnum; else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; @@ -627,10 +629,10 @@ else if (strcmp(name,"VzMesh")==0) return VzMeshEnum; else if (strcmp(name,"VzSSA")==0) return VzSSAEnum; - else if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum; - else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum; + if (strcmp(name,"Watercolumn")==0) return WatercolumnEnum; + else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum; + else if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum; else if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum; else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum; @@ -750,10 +752,10 @@ else if (strcmp(name,"FSpressure")==0) return FSpressureEnum; else if (strcmp(name,"FSSolver")==0) return FSSolverEnum; - else if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum; - else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; + if (strcmp(name,"FSvelocity")==0) return FSvelocityEnum; + else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum; + else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; @@ -873,10 +875,10 @@ else if (strcmp(name,"MaxAbsVz")==0) return MaxAbsVzEnum; else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum; - else if (strcmp(name,"MaxVel")==0) return MaxVelEnum; - else if (strcmp(name,"MaxVx")==0) return MaxVxEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"MaxVy")==0) return MaxVyEnum; + if (strcmp(name,"MaxVel")==0) return MaxVelEnum; + else if (strcmp(name,"MaxVx")==0) return MaxVxEnum; + else if (strcmp(name,"MaxVy")==0) return MaxVyEnum; else if (strcmp(name,"MaxVz")==0) return MaxVzEnum; else if (strcmp(name,"Melange")==0) return MelangeEnum; @@ -996,10 +998,10 @@ else if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum; else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum; - else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum; - else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum; else stage=9; } if(stage==9){ - if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum; + if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum; + else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum; + else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum; else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum; else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum; @@ -1119,10 +1121,10 @@ else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum; else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; - else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; - else if (strcmp(name,"TransientParam")==0) return TransientParamEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; + if (strcmp(name,"TransientInput")==0) return TransientInputEnum; + else if (strcmp(name,"TransientParam")==0) return TransientParamEnum; + else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; else if (strcmp(name,"Tria")==0) return TriaEnum; else if (strcmp(name,"TriaInput")==0) return TriaInputEnum; Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 22897) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp (revision 22898) @@ -1,5 +1,5 @@ /*!\file: solutionsequence_hydro_nonlinear.cpp - * \brief: core of the hydro solution - */ + * \brief: core of the hydro solution + */ #include "../toolkits/toolkits.h" @@ -12,16 +12,16 @@ void solutionsequence_hydro_nonlinear(FemModel* femmodel){ /*solution : */ - Vector* ug_sed=NULL; - Vector* uf_sed=NULL; - Vector* uf_sed_sub_iter=NULL; + Vector* ug_sed=NULL; + Vector* uf_sed=NULL; + Vector* uf_sed_sub_iter=NULL; Vector* ug_sed_main_iter=NULL; - Vector* ug_epl=NULL; + Vector* ug_epl=NULL; Vector* uf_epl=NULL; - Vector* uf_epl_sub_iter=NULL; + Vector* uf_epl_sub_iter=NULL; Vector* ug_epl_sub_iter=NULL; Vector* ug_epl_main_iter=NULL; - Vector* ys=NULL; + Vector* ys=NULL; Vector* dug=NULL; Vector* duf=NULL; @@ -35,5 +35,5 @@ HydrologyDCInefficientAnalysis* inefanalysis = NULL; HydrologyDCEfficientAnalysis* effanalysis = NULL; - + bool sedconverged,eplconverged,hydroconverged; bool isefficientlayer; @@ -58,6 +58,7 @@ /*Retrieve inputs as the initial state for the non linear iteration*/ - GetSolutionFromInputsx(&ug_sed,femmodel); + GetSolutionFromInputsx(&ug_sed,femmodel); Reducevectorgtofx(&uf_sed, ug_sed, femmodel->nodes,femmodel->parameters); + /*Initialize the thawed element mask*/ if(isefficientlayer) { inefanalysis = new HydrologyDCInefficientAnalysis(); @@ -87,5 +88,5 @@ InputUpdateFromConstantx(femmodel,false,ConvergedEnum); femmodel->UpdateConstraintsx(); - + /*Reset constraint on the ZigZag Lock*/ ResetConstraintsx(femmodel); @@ -126,8 +127,8 @@ if(sedconverged)break; } - + /*}}}*//*End of the sediment penalization loop*/ sedconverged=false; - + /*Checking convegence on the value of the sediment head*/ duf=uf_sed_sub_iter->Duplicate();_assert_(duf); @@ -184,5 +185,5 @@ inefanalysis->ElementizeEplMask(femmodel); /*}}}*/ - + if(VerboseSolution()) _printf0_("Building EPL Matrix...\n"); SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel); @@ -197,9 +198,9 @@ uf_epl_sub_iter=uf_epl->Duplicate(); uf_epl->Copy(uf_epl_sub_iter); - delete ug_epl; + delete ug_epl; Mergesolutionfromftogx(&ug_epl,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys; InputUpdateFromSolutionx(femmodel,ug_epl); ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); - + dug=ug_epl_sub_iter->Duplicate();_assert_(dug); ug_epl_sub_iter->Copy(dug); @@ -218,5 +219,5 @@ /* if(ThickCountDuplicate(); _assert_(dug); - ug_sed_main_iter->Copy(dug); + ug_sed_main_iter->Copy(dug); dug->AYPX(ug_sed,-1.0); - ndu_sed=dug->Norm(NORM_TWO); + ndu_sed=dug->Norm(NORM_TWO); delete dug; nu_sed=ug_sed_main_iter->Norm(NORM_TWO); @@ -243,8 +244,8 @@ if (xIsNan(ndu_sed) || xIsNan(nu_sed)) _error_("Sed convergence criterion is NaN!"); if (ndu_sed==0.0 && nu_sed==0.0) nu_sed=1.0e-6; /*Hacking the case where the Sediment is used but empty*/ - dug=ug_epl_main_iter->Duplicate();_assert_(dug); - ug_epl_main_iter->Copy(dug); + dug=ug_epl_main_iter->Duplicate();_assert_(dug); + ug_epl_main_iter->Copy(dug); dug->AYPX(ug_epl,-1.0); - ndu_epl=dug->Norm(NORM_TWO); + ndu_epl=dug->Norm(NORM_TWO); delete dug; nu_epl=ug_epl_main_iter->Norm(NORM_TWO); @@ -257,5 +258,5 @@ hydroconverged=true; } - else{ + else{ if(VerboseConvergence()) _printf0_(setw(50) << left << " Sediment Convergence criterion:" << ndu_sed/nu_sed*100 << "%, aiming lower than " << eps_hyd*100 << " %\n"); if(VerboseConvergence()) _printf0_(setw(50) << left << " EPL Convergence criterion:" << ndu_epl/nu_epl*100 << "%, aiming lower than " << eps_hyd*100 << " %\n");