Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 24090) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 24091) @@ -624,19 +624,14 @@ void DamageEvolutionAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - int domaintype; - IssmDouble max_damage; - int *doflist = NULL; - - element->FindParam(&domaintype,DomainTypeEnum); - /*Fetch number of nodes and dof for this finite element*/ int numnodes = element->GetNumberOfNodes(); /*Fetch dof list and allocate solution vector*/ + int *doflist = NULL; element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum); IssmDouble* newdamage = xNew(numnodes); /*Get user-supplied max_damage: */ - element->FindParam(&max_damage,DamageMaxDamageEnum); + IssmDouble max_damage = element->FindParam(DamageMaxDamageEnum); /*Use the dof list to index into the solution vector: */ @@ -652,4 +647,6 @@ /*Get all inputs and parameters*/ + int domaintype; + element->FindParam(&domaintype,DomainTypeEnum); if(domaintype==Domain2DhorizontalEnum){ element->AddInput(DamageDbarEnum,newdamage,element->GetElementType()); Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 24090) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 24091) @@ -5,4 +5,6 @@ #include "../modules/modules.h" +#define FINITEELEMENT P1Enum + /*Model processing*/ void MasstransportAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ @@ -14,5 +16,5 @@ /*Do not add constraints in DG, they are weakly imposed*/ if(stabilization!=3){ - IoModelToConstraintsx(constraints,iomodel,"md.masstransport.spcthickness",MasstransportAnalysisEnum,P1Enum); + IoModelToConstraintsx(constraints,iomodel,"md.masstransport.spcthickness",MasstransportAnalysisEnum,FINITEELEMENT); } @@ -25,5 +27,4 @@ /*Intermediaries*/ - int element; int penpair_ids[2]; int count=0; @@ -35,5 +36,5 @@ /*Loads only in DG*/ - if (stabilization==3){ + if(stabilization==3){ /*Get faces and elements*/ @@ -45,10 +46,8 @@ /*Get left and right elements*/ - element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2] + int element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2] /*Now, if this element is not in the partition, pass: */ if(!iomodel->my_elements[element]) continue; - - /* Add load */ loads->AddObject(new Numericalflux(i+1,i,i,iomodel)); } @@ -102,5 +101,5 @@ if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); if(stabilization!=3){ - ::CreateNodes(nodes,iomodel,MasstransportAnalysisEnum,P1Enum,isamr); + ::CreateNodes(nodes,iomodel,MasstransportAnalysisEnum,FINITEELEMENT,isamr); } else{ @@ -130,7 +129,8 @@ /*Finite element type*/ - finiteelement = P1Enum; + finiteelement = FINITEELEMENT; if(stabilization==3){ finiteelement = P1DGEnum; + //finiteelement = P0DGEnum; } @@ -260,4 +260,5 @@ Ke = CreateKMatrixCG(basalelement); break; + case P0DGEnum: case P1DGEnum: Ke = CreateKMatrixDG(basalelement); @@ -549,4 +550,5 @@ pe = CreatePVectorCG(basalelement); break; + case P0DGEnum: case P1DGEnum: pe = CreatePVectorDG(basalelement); @@ -706,9 +708,9 @@ element->GetVerticesCoordinates(&xyz_list); element->FindParam(&dt,TimesteppingTimeStepEnum); - Input* gmb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input); - Input* fmb_input = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input); - Input* ms_input = element->GetInput(SmbMassBalanceEnum); _assert_(ms_input); - Input* gllevelset_input = element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input); - Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* gmb_input = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input); + Input* fmb_input = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input); + Input* ms_input = element->GetInput(SmbMassBalanceEnum); _assert_(ms_input); + Input* gllevelset_input = element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input); + Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); /* Start looping on the number of gaussian points: */ @@ -803,59 +805,68 @@ void MasstransportAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - int i,hydroadjustment,domaintype; - int* doflist=NULL; - IssmDouble rho_ice,rho_water,minthickness; - Element* basalelement=NULL; - bool isgroundingline=0; - - element->FindParam(&domaintype,DomainTypeEnum); - if(domaintype!=Domain2DhorizontalEnum){ - if(!element->IsOnBase()) return; - basalelement=element->SpawnBasalElement(); - } - else{ - basalelement = element; - } - - /*Fetch number of nodes for this finite element*/ - int numnodes = basalelement->GetNumberOfNodes(); + /*Only update if on base*/ + if(!element->IsOnBase()) return; + + /*Get basal element*/ + int domaintype; element->FindParam(&domaintype,DomainTypeEnum); + Element* basalelement=element; + if(domaintype!=Domain2DhorizontalEnum) basalelement = element->SpawnBasalElement(); + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = basalelement->GetNumberOfNodes(); + int numvertices = basalelement->GetNumberOfVertices(); + + /*Keep old thickness for later*/ + IssmDouble* oldthickness = xNew(numvertices); + basalelement->GetInputListOnNodes(&oldthickness[0],ThicknessEnum); /*Fetch dof list and allocate solution vector*/ - basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum); - IssmDouble* newthickness = xNew(numnodes); - IssmDouble* cumdeltathickness = xNew(numnodes); - IssmDouble* deltathickness = xNew(numnodes); - IssmDouble* newbase = xNew(numnodes); - IssmDouble* bed = xNew(numnodes); - IssmDouble* newsurface = xNew(numnodes); - IssmDouble* oldthickness = xNew(numnodes); - IssmDouble* oldbase = xNew(numnodes); - IssmDouble* oldsurface = xNew(numnodes); - IssmDouble* phi = xNew(numnodes); - IssmDouble* sealevel = xNew(numnodes); + int *doflist = NULL; + element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum); + IssmDouble* newthickness = xNew(numnodes); /*Use the dof list to index into the solution vector: */ - basalelement->FindParam(&minthickness,MasstransportMinThicknessEnum); - for(i=0;iFindParam(MasstransportMinThicknessEnum); + for(int i=0;i(newthickness[i])) _error_("NaN found in solution vector"); if(xIsInf(newthickness[i])) _error_("Inf found in solution vector"); - /*Constrain thickness to be at least 1m*/ - if(newthickness[i]AddBasalInput(ThicknessEnum,newthickness,element->GetElementType()); + + /*Free ressources:*/ + xDelete(newthickness); + xDelete(doflist); + + /*Now, we need to do some "processing"*/ + newthickness = xNew(numvertices); + IssmDouble* cumdeltathickness = xNew(numvertices); + IssmDouble* deltathickness = xNew(numvertices); + IssmDouble* newbase = xNew(numvertices); + IssmDouble* bed = xNew(numvertices); + IssmDouble* newsurface = xNew(numvertices); + IssmDouble* oldbase = xNew(numvertices); + IssmDouble* oldsurface = xNew(numvertices); + IssmDouble* phi = xNew(numvertices); + IssmDouble* sealevel = xNew(numvertices); /*Get previous base, thickness, surfac and current sealevel and bed:*/ - basalelement->GetInputListOnNodes(&oldbase[0],BaseEnum); - basalelement->GetInputListOnNodes(&oldsurface[0],SurfaceEnum); - basalelement->GetInputListOnNodes(&oldthickness[0],ThicknessEnum); - basalelement->GetInputListOnNodes(&phi[0],MaskGroundediceLevelsetEnum); - basalelement->GetInputListOnNodes(&sealevel[0],SealevelEnum); - basalelement->GetInputListOnNodes(&cumdeltathickness[0],SealevelriseCumDeltathicknessEnum); - + basalelement->GetInputListOnVertices(&newthickness[0],ThicknessEnum); + basalelement->GetInputListOnVertices(&oldbase[0],BaseEnum); + basalelement->GetInputListOnVertices(&oldsurface[0],SurfaceEnum); + basalelement->GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum); + basalelement->GetInputListOnVertices(&sealevel[0],SealevelEnum); + basalelement->GetInputListOnVertices(&cumdeltathickness[0],SealevelriseCumDeltathicknessEnum); + + /*Do we do grounding line migration?*/ + bool isgroundingline; element->FindParam(&isgroundingline,TransientIsgroundinglineEnum); - if(isgroundingline) basalelement->GetInputListOnNodes(&bed[0],BedEnum); + if(isgroundingline) basalelement->GetInputListOnVertices(&bed[0],BedEnum); /*What is the delta thickness forcing the sea-level rise core: cumulated over time, hence the +=:*/ - for(i=0;iFindParam(&hydroadjustment,MasstransportHydrostaticAdjustmentEnum); - rho_ice = basalelement->FindParam(MaterialsRhoIceEnum); - rho_water = basalelement->FindParam(MaterialsRhoSeawaterEnum); - - for(i=0;iFindParam(MaterialsRhoIceEnum); + IssmDouble rho_water = basalelement->FindParam(MaterialsRhoSeawaterEnum); + + for(int i=0;i0.){ //this is an ice sheet: just add thickness to base. /*Update! actually, the bed has moved too!:*/ if(isgroundingline){ newsurface[i] = bed[i]+newthickness[i]; //surface = bed + newthickness - newbase[i] = bed[i]; //new base at new bed + newbase[i] = bed[i]; //new base at new bed } else{ @@ -881,10 +893,10 @@ else{ //this is an ice shelf: hydrostatic equilibrium*/ if(hydroadjustment==AbsoluteEnum){ - newsurface[i] = newthickness[i]*(1-rho_ice/rho_water)+sealevel[i]; - newbase[i] = newthickness[i]*(-rho_ice/rho_water)+sealevel[i]; + newsurface[i] = newthickness[i]*(1.-rho_ice/rho_water)+sealevel[i]; + newbase[i] = newthickness[i]*(-rho_ice/rho_water)+sealevel[i]; } else if(hydroadjustment==IncrementalEnum){ newsurface[i] = oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i])+sealevel[i]; //surface = oldsurface + (1-di) * dH - newbase[i] = oldbase[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i])+sealevel[i]; //base = oldbed + di * dH + newbase[i] = oldbase[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i])+sealevel[i]; //base = oldbed + di * dH } else _error_("Hydrostatic adjustment " << hydroadjustment << " (" << EnumToStringx(hydroadjustment) << ") not supported yet"); @@ -893,5 +905,4 @@ /*Add input to the element: */ - element->AddBasalInput(ThicknessEnum,newthickness,P1Enum); element->AddBasalInput(SurfaceEnum,newsurface,P1Enum); element->AddBasalInput(BaseEnum,newbase,P1Enum); @@ -911,5 +922,4 @@ xDelete(sealevel); xDelete(bed); - xDelete(doflist); if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24090) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24091) @@ -20,5 +20,5 @@ /* SEMIC prototype {{{*/ extern "C" void run_semic_(double *sf_in, double *rf_in, double *swd_in, double *lwd_in, double *wind_in, double *sp_in, double *rhoa_in, - double *qq_in, double *tt_in, double *tsurf_out, double *smb_out, double *saccu_out, double *smelt_out); + double *qq_in, double *tt_in, double *tsurf_out, double *smb_out, double *saccu_out, double *smelt_out); #endif // _HAVE_SEMIC_ @@ -161,21 +161,21 @@ /* Fetch number of nodes and allocate output*/ - int numnodes = this->GetNumberOfNodes(); - IssmDouble* newD = xNew(numnodes); + int numnodes = this->GetNumberOfNodes(); + IssmDouble* newD = xNew(numnodes); /* Retrieve domain-dependent inputs */ Input* n_input=this->GetInput(MaterialsRheologyNEnum); _assert_(n_input); - Input* damage_input = NULL; - Input* B_input = NULL; + Input* damage_input = NULL; + Input* B_input = NULL; int domaintype; - parameters->FindParam(&domaintype,DomainTypeEnum); + parameters->FindParam(&domaintype,DomainTypeEnum); if(domaintype==Domain2DhorizontalEnum){ - damage_input = this->GetInput(DamageDbarEnum); _assert_(damage_input); - B_input=this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input); - } - else{ - damage_input = this->GetInput(DamageDEnum); _assert_(damage_input); - B_input=this->GetInput(MaterialsRheologyBEnum); _assert_(B_input); - } + damage_input = this->GetInput(DamageDbarEnum); _assert_(damage_input); + B_input=this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input); + } + else{ + damage_input = this->GetInput(DamageDEnum); _assert_(damage_input); + B_input=this->GetInput(MaterialsRheologyBEnum); _assert_(B_input); + } /* Start looping on the number of nodes: */ @@ -198,6 +198,6 @@ B_input->GetInputValue(&B,gauss); - n_input->GetInputValue(&n,gauss); - damage_input->GetInputValue(&D,gauss); + n_input->GetInputValue(&n,gauss); + damage_input->GetInputValue(&D,gauss); /* Compute threshold strain rate from threshold stress */ @@ -268,5 +268,5 @@ if(dim==2){ - /* epsilon=[exx,eyy,exy];*/ + /* epsilon=[exx,eyy,exy];*/ eps_xx[iv]=epsilon[0]; eps_yy[iv]=epsilon[1]; @@ -705,21 +705,21 @@ } switch(this->ObjectEnum()){ - case TriaEnum: - this->inputs->AddInput(new TriaInput(SmbMassBalanceEnum,&smb[0],P1Enum)); - this->inputs->AddInput(new TriaInput(SmbRunoffEnum,&runoff[0],P1Enum)); - break; - case PentaEnum: - this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&smb[0],P1Enum)); - this->inputs->AddInput(new PentaInput(SmbRunoffEnum,&runoff[0],P1Enum)); - this->InputExtrude(SmbMassBalanceEnum,-1); - this->InputExtrude(SmbRunoffEnum,-1); - break; - case TetraEnum: - this->inputs->AddInput(new TetraInput(SmbMassBalanceEnum,&smb[0],P1Enum)); - this->inputs->AddInput(new TetraInput(SmbRunoffEnum,&runoff[0],P1Enum)); - this->InputExtrude(SmbMassBalanceEnum,-1); - this->InputExtrude(SmbRunoffEnum,-1); - break; - default: _error_("Not implemented yet"); + case TriaEnum: + this->inputs->AddInput(new TriaInput(SmbMassBalanceEnum,&smb[0],P1Enum)); + this->inputs->AddInput(new TriaInput(SmbRunoffEnum,&runoff[0],P1Enum)); + break; + case PentaEnum: + this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&smb[0],P1Enum)); + this->inputs->AddInput(new PentaInput(SmbRunoffEnum,&runoff[0],P1Enum)); + this->InputExtrude(SmbMassBalanceEnum,-1); + this->InputExtrude(SmbRunoffEnum,-1); + break; + case TetraEnum: + this->inputs->AddInput(new TetraInput(SmbMassBalanceEnum,&smb[0],P1Enum)); + this->inputs->AddInput(new TetraInput(SmbRunoffEnum,&runoff[0],P1Enum)); + this->InputExtrude(SmbMassBalanceEnum,-1); + this->InputExtrude(SmbRunoffEnum,-1); + break; + default: _error_("Not implemented yet"); } /* this->AddInput(SmbMassBalanceEnum,smb,P1Enum); */ @@ -858,14 +858,14 @@ IssmDouble eps_eff; - if(dim==2){ - /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/ - this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); - eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]); - } - else{ - /* eps_eff^2 = 1/2 exx^2*/ - this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input); - eps_eff = sqrt(epsilon1d*epsilon1d/2.); - } + if(dim==2){ + /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/ + this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input); + eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]); + } + else{ + /* eps_eff^2 = 1/2 exx^2*/ + this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input); + eps_eff = sqrt(epsilon1d*epsilon1d/2.); + } /*Get viscosity*/ @@ -1319,7 +1319,7 @@ }/*}}}*/ Node* Element::GetNode(int nodeindex){/*{{{*/ - _assert_(nodeindex>=0); - _assert_(nodeindexGetNumberOfNodes(this->element_type)); - return this->nodes[nodeindex]; + _assert_(nodeindex>=0); + _assert_(nodeindexGetNumberOfNodes(this->element_type)); + return this->nodes[nodeindex]; }/*}}}*/ int Element::GetNodeIndex(Node* node){/*{{{*/ @@ -1511,22 +1511,22 @@ switch(type){ - case VertexPIdEnum: - doflist = xNew(NUM_VERTICES); - values = xNew(NUM_VERTICES); - /*Fill in values*/ - this->GetVerticesPidList(doflist); - this->GetInputListOnVerticesAtTime(values,input_enum,time); - vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL); - break; - case VertexSIdEnum: - doflist = xNew(NUM_VERTICES); - values = xNew(NUM_VERTICES); - /*Fill in values*/ - this->GetVerticesSidList(doflist); - this->GetInputListOnVerticesAtTime(values,input_enum,time); - vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL); - break; - default: - _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); + case VertexPIdEnum: + doflist = xNew(NUM_VERTICES); + values = xNew(NUM_VERTICES); + /*Fill in values*/ + this->GetVerticesPidList(doflist); + this->GetInputListOnVerticesAtTime(values,input_enum,time); + vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL); + break; + case VertexSIdEnum: + doflist = xNew(NUM_VERTICES); + values = xNew(NUM_VERTICES); + /*Fill in values*/ + this->GetVerticesSidList(doflist); + this->GetInputListOnVerticesAtTime(values,input_enum,time); + vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL); + break; + default: + _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); } @@ -1741,119 +1741,119 @@ void Element::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){/*{{{*/ - /*Intermediaries*/ - int i,t; - - /*Branch on type of vector: nodal or elementary: */ - if(vector_type==1){ //nodal vector - - const int NUM_VERTICES = this->GetNumberOfVertices(); - - int *vertexids = xNew(NUM_VERTICES); - IssmDouble *values = xNew(NUM_VERTICES); - - /*Recover vertices ids needed to initialize inputs*/ - _assert_(iomodel->elements); - for(i=0;i(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab - } - - /*Are we in transient or static? */ - if(M==1){ - values[0]=vector[0]; - this->AddInput(vector_enum,values,P0Enum); - } - else if(M==iomodel->numberofvertices){ - for(i=0;iAddInput(vector_enum,values,P1Enum); - } - else if(M==iomodel->numberofvertices+1){ - /*create transient input: */ - IssmDouble* times = xNew(N); - for(t=0;tObjectEnum()){ - case TriaEnum: transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break; - case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,values,P1Enum)); break; - case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,values,P1Enum)); break; - default: _error_("Not implemented yet"); - } - } - this->inputs->AddInput(transientinput); - xDelete(times); - } - else if(M==iomodel->numberofelements){ - - /*This is a Patch!*/ - xDelete(values); - values = xNew(N); - for(int j=0;jSid()*N+j]; - - if (N==this->GetNumberOfNodes(P1Enum) ) this->AddInput(vector_enum,values,P1Enum); - else if(N==this->GetNumberOfNodes(P0Enum) ) this->AddInput(vector_enum,values,P0Enum); - else if(N==this->GetNumberOfNodes(P1xP2Enum)) this->AddInput(vector_enum,values,P1xP2Enum); - else if(N==this->GetNumberOfNodes(P1xP3Enum)) this->AddInput(vector_enum,values,P1xP3Enum); - else _error_("Patch interpolation not supported yet"); - - } - else{ - _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); - } - - xDelete(values); - xDelete(vertexids); - } - else if(vector_type==2){ //element vector - - IssmDouble value; - - /*Are we in transient or static? */ - if(M==iomodel->numberofelements){ - if (code==5){ //boolean - this->inputs->AddInput(new BoolInput(vector_enum,reCast(vector[this->Sid()]))); - } - else if (code==6){ //integer - this->inputs->AddInput(new IntInput(vector_enum,reCast(vector[this->Sid()]))); - } - else if (code==7){ //IssmDouble - this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->Sid()])); - } - else _error_("could not recognize nature of vector from code " << code); - } - else if(M==iomodel->numberofelements+1){ - /*create transient input: */ - IssmDouble* times = xNew(N); - for(t=0;tSid()+t]; - switch(this->ObjectEnum()){ - case TriaEnum: transientinput->AddTimeInput(new TriaInput( vector_enum,&value,P0Enum)); break; - case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,&value,P0Enum)); break; - case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,&value,P0Enum)); break; - default: _error_("Not implemented yet"); - } - } - this->inputs->AddInput(transientinput); - xDelete(times); - } - else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); - } - else if(vector_type==3){ //element vector - - /*For right now we are static */ - if(M==iomodel->numberofelements){ - /*create transient input: */ - IssmDouble* layers = xNewZeroInit(N);; - for(t=0;tSid()+t]; - DoubleArrayInput* arrayinput=new DoubleArrayInput(vector_enum,layers,N); - this->inputs->AddInput(arrayinput); - xDelete(layers); - } - else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); - } - else _error_("Cannot add input for vector type " << vector_type << " (not supported)"); + /*Intermediaries*/ + int i,t; + + /*Branch on type of vector: nodal or elementary: */ + if(vector_type==1){ //nodal vector + + const int NUM_VERTICES = this->GetNumberOfVertices(); + + int *vertexids = xNew(NUM_VERTICES); + IssmDouble *values = xNew(NUM_VERTICES); + + /*Recover vertices ids needed to initialize inputs*/ + _assert_(iomodel->elements); + for(i=0;i(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab + } + + /*Are we in transient or static? */ + if(M==1){ + values[0]=vector[0]; + this->AddInput(vector_enum,values,P0Enum); + } + else if(M==iomodel->numberofvertices){ + for(i=0;iAddInput(vector_enum,values,P1Enum); + } + else if(M==iomodel->numberofvertices+1){ + /*create transient input: */ + IssmDouble* times = xNew(N); + for(t=0;tObjectEnum()){ + case TriaEnum: transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break; + case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,values,P1Enum)); break; + case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,values,P1Enum)); break; + default: _error_("Not implemented yet"); + } + } + this->inputs->AddInput(transientinput); + xDelete(times); + } + else if(M==iomodel->numberofelements){ + + /*This is a Patch!*/ + xDelete(values); + values = xNew(N); + for(int j=0;jSid()*N+j]; + + if (N==this->GetNumberOfNodes(P1Enum) ) this->AddInput(vector_enum,values,P1Enum); + else if(N==this->GetNumberOfNodes(P0Enum) ) this->AddInput(vector_enum,values,P0Enum); + else if(N==this->GetNumberOfNodes(P1xP2Enum)) this->AddInput(vector_enum,values,P1xP2Enum); + else if(N==this->GetNumberOfNodes(P1xP3Enum)) this->AddInput(vector_enum,values,P1xP3Enum); + else _error_("Patch interpolation not supported yet"); + + } + else{ + _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); + } + + xDelete(values); + xDelete(vertexids); + } + else if(vector_type==2){ //element vector + + IssmDouble value; + + /*Are we in transient or static? */ + if(M==iomodel->numberofelements){ + if (code==5){ //boolean + this->inputs->AddInput(new BoolInput(vector_enum,reCast(vector[this->Sid()]))); + } + else if (code==6){ //integer + this->inputs->AddInput(new IntInput(vector_enum,reCast(vector[this->Sid()]))); + } + else if (code==7){ //IssmDouble + this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->Sid()])); + } + else _error_("could not recognize nature of vector from code " << code); + } + else if(M==iomodel->numberofelements+1){ + /*create transient input: */ + IssmDouble* times = xNew(N); + for(t=0;tSid()+t]; + switch(this->ObjectEnum()){ + case TriaEnum: transientinput->AddTimeInput(new TriaInput( vector_enum,&value,P0Enum)); break; + case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,&value,P0Enum)); break; + case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,&value,P0Enum)); break; + default: _error_("Not implemented yet"); + } + } + this->inputs->AddInput(transientinput); + xDelete(times); + } + else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); + } + else if(vector_type==3){ //element vector + + /*For right now we are static */ + if(M==iomodel->numberofelements){ + /*create transient input: */ + IssmDouble* layers = xNewZeroInit(N);; + for(t=0;tSid()+t]; + DoubleArrayInput* arrayinput=new DoubleArrayInput(vector_enum,layers,N); + this->inputs->AddInput(arrayinput); + xDelete(layers); + } + else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); + } + else _error_("Cannot add input for vector type " << vector_type << " (not supported)"); } /*}}}*/ @@ -1894,41 +1894,41 @@ else if(M==iomodel->numberofvertices+1){ - /*create transient input: */ - IssmDouble* times = xNew(N); - for(int t=0;tObjectEnum()){ - case TriaEnum: - values_input->AddTimeInput(new TriaInput(input_enum,values,P1Enum)); - mins_input->AddTimeInput(new TriaInput(ControlInputMinsEnum,values_min,P1Enum)); - maxs_input->AddTimeInput(new TriaInput(ControlInputMaxsEnum,values_max,P1Enum)); - break; - case PentaEnum: - values_input->AddTimeInput(new PentaInput(input_enum,values,P1Enum)); - mins_input->AddTimeInput(new PentaInput(ControlInputMinsEnum,values_min,P1Enum)); - maxs_input->AddTimeInput(new PentaInput(ControlInputMaxsEnum,values_max,P1Enum)); - break; - case TetraEnum: - values_input->AddTimeInput(new TetraInput(input_enum,values,P1Enum)); - mins_input->AddTimeInput(new TetraInput(ControlInputMinsEnum,values_min,P1Enum)); - maxs_input->AddTimeInput(new TetraInput(ControlInputMaxsEnum,values_max,P1Enum)); - break; - default: _error_("Not implemented yet"); - } - } - this->inputs->AddInput(new ControlInput(input_enum,TransientInputEnum,values_input,mins_input,maxs_input,grad_input,P1Enum,id)); - xDelete(times); - } - else _error_("not currently supported type of M and N attempted"); + /*create transient input: */ + IssmDouble* times = xNew(N); + for(int t=0;tObjectEnum()){ + case TriaEnum: + values_input->AddTimeInput(new TriaInput(input_enum,values,P1Enum)); + mins_input->AddTimeInput(new TriaInput(ControlInputMinsEnum,values_min,P1Enum)); + maxs_input->AddTimeInput(new TriaInput(ControlInputMaxsEnum,values_max,P1Enum)); + break; + case PentaEnum: + values_input->AddTimeInput(new PentaInput(input_enum,values,P1Enum)); + mins_input->AddTimeInput(new PentaInput(ControlInputMinsEnum,values_min,P1Enum)); + maxs_input->AddTimeInput(new PentaInput(ControlInputMaxsEnum,values_max,P1Enum)); + break; + case TetraEnum: + values_input->AddTimeInput(new TetraInput(input_enum,values,P1Enum)); + mins_input->AddTimeInput(new TetraInput(ControlInputMinsEnum,values_min,P1Enum)); + maxs_input->AddTimeInput(new TetraInput(ControlInputMaxsEnum,values_max,P1Enum)); + break; + default: _error_("Not implemented yet"); + } + } + this->inputs->AddInput(new ControlInput(input_enum,TransientInputEnum,values_input,mins_input,maxs_input,grad_input,P1Enum,id)); + xDelete(times); + } + else _error_("not currently supported type of M and N attempted"); /*clean up*/ @@ -1938,176 +1938,176 @@ /*}}}*/ void Element::DatasetInputAdd(int enum_type,IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int input_enum,int code,int input_id){/*{{{*/ - /*enum_type: the name of the DatasetInput (eg Outputdefinition1) - * vector: information being stored (eg observations) - * vector_type: is if by element or by vertex - * input_enum: is the name of the vector being stored - * code: what type of data is in the vector (booleans, ints, doubles) - */ - - /*Intermediaries*/ - int i,t; - DatasetInput* datasetinput = NULL; - - /*Get input if it already exists*/ - Input* tempinput = GetInput(enum_type); - if(tempinput){ - /*Cast it to a Datasetinput*/ - if(tempinput->ObjectEnum()!=DatasetInputEnum) _error_("don't know what to do"); - datasetinput = (DatasetInput*)tempinput; - } - else{ - datasetinput=new DatasetInput(enum_type); - this->inputs->AddInput(datasetinput); - } - - /*Branch on type of vector: nodal or elementary: */ - if(vector_type==1){ //nodal vector - - const int NUM_VERTICES = this->GetNumberOfVertices(); - - int *vertexids = xNew(NUM_VERTICES); - IssmDouble *values = xNew(NUM_VERTICES); - - /*Recover vertices ids needed to initialize inputs*/ - _assert_(iomodel->elements); - for(i=0;i(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab - } - - /*Are we in transient or static? */ - if(M==1){ - values[0]=vector[0]; + /*enum_type: the name of the DatasetInput (eg Outputdefinition1) + * vector: information being stored (eg observations) + * vector_type: is if by element or by vertex + * input_enum: is the name of the vector being stored + * code: what type of data is in the vector (booleans, ints, doubles) + */ + + /*Intermediaries*/ + int i,t; + DatasetInput* datasetinput = NULL; + + /*Get input if it already exists*/ + Input* tempinput = GetInput(enum_type); + if(tempinput){ + /*Cast it to a Datasetinput*/ + if(tempinput->ObjectEnum()!=DatasetInputEnum) _error_("don't know what to do"); + datasetinput = (DatasetInput*)tempinput; + } + else{ + datasetinput=new DatasetInput(enum_type); + this->inputs->AddInput(datasetinput); + } + + /*Branch on type of vector: nodal or elementary: */ + if(vector_type==1){ //nodal vector + + const int NUM_VERTICES = this->GetNumberOfVertices(); + + int *vertexids = xNew(NUM_VERTICES); + IssmDouble *values = xNew(NUM_VERTICES); + + /*Recover vertices ids needed to initialize inputs*/ + _assert_(iomodel->elements); + for(i=0;i(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab + } + + /*Are we in transient or static? */ + if(M==1){ + values[0]=vector[0]; + switch(this->ObjectEnum()){ + case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P0Enum),input_id); break; + case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P0Enum),input_id); break; + case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P0Enum),input_id); break; + default: _error_("Not implemented yet"); + } + } + else if(M==iomodel->numberofvertices){ + for(i=0;iObjectEnum()){ + case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P1Enum),input_id); break; + case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1Enum),input_id); break; + case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1Enum),input_id); break; + default: _error_("Not implemented yet"); + } } + else if(M==iomodel->numberofvertices+1){ + /*create transient input: */ + IssmDouble* times = xNew(N); + for(t=0;tObjectEnum()){ - case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P0Enum),input_id); break; - case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P0Enum),input_id); break; - case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P0Enum),input_id); break; - default: _error_("Not implemented yet"); + case TriaEnum: transientinput->AddTimeInput(new TriaInput( input_enum,values,P1Enum)); break; + case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,values,P1Enum)); break; + case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,values,P1Enum)); break; + default: _error_("Not implemented yet"); } - } - else if(M==iomodel->numberofvertices){ - for(i=0;iAddInput(transientinput,input_id); + xDelete(times); + } + else if(M==iomodel->numberofelements){ + + /*This is a Patch!*/ + xDelete(values); + values = xNew(N); + for(int j=0;jSid()*N+j]; + + if (N==this->GetNumberOfNodes(P1Enum) ){ switch(this->ObjectEnum()){ - case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P1Enum),input_id); break; - case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1Enum),input_id); break; - case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1Enum),input_id); break; - default: _error_("Not implemented yet"); - } } - else if(M==iomodel->numberofvertices+1){ - /*create transient input: */ - IssmDouble* times = xNew(N); - for(t=0;tObjectEnum()){ - case TriaEnum: transientinput->AddTimeInput(new TriaInput( input_enum,values,P1Enum)); break; - case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,values,P1Enum)); break; - case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,values,P1Enum)); break; - default: _error_("Not implemented yet"); - } - } - datasetinput->AddInput(transientinput,input_id); - xDelete(times); - } - else if(M==iomodel->numberofelements){ - - /*This is a Patch!*/ - xDelete(values); - values = xNew(N); - for(int j=0;jSid()*N+j]; - - if (N==this->GetNumberOfNodes(P1Enum) ){ - switch(this->ObjectEnum()){ - case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P1Enum),input_id); break; - case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1Enum),input_id); break; - case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1Enum),input_id); break; - default: _error_("Not implemented yet"); - } - } - else if(N==this->GetNumberOfNodes(P0Enum) ){ - switch(this->ObjectEnum()){ - case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P0Enum),input_id); break; - case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P0Enum),input_id); break; - case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P0Enum),input_id); break; - default: _error_("Not implemented yet"); - } - } - else if(N==this->GetNumberOfNodes(P1xP2Enum)){ - switch(this->ObjectEnum()){ - case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P1xP2Enum),input_id); break; - case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1xP2Enum),input_id); break; - case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1xP2Enum),input_id); break; - default: _error_("Not implemented yet"); - } - } - else if(N==this->GetNumberOfNodes(P1xP3Enum)) { - switch(this->ObjectEnum()){ - case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P1xP3Enum),input_id); break; - case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1xP3Enum),input_id); break; - case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1xP3Enum),input_id); break; - default: _error_("Not implemented yet"); - } - } - else _error_("Patch interpolation not supported yet"); - - } - else{ - _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long"); - } - - xDelete(values); - xDelete(vertexids); - } - else if(vector_type==2){ //element vector - - IssmDouble value; - - /*Are we in transient or static? */ - if(M==iomodel->numberofelements){ - if (code==5){ //boolean - datasetinput->AddInput(new BoolInput(input_enum,reCast(vector[this->Sid()])),input_id); - } - else if (code==6){ //integer - datasetinput->AddInput(new IntInput(input_enum,reCast(vector[this->Sid()])),input_id); - } - else if (code==7){ //IssmDouble - datasetinput->AddInput(new DoubleInput(input_enum,vector[this->Sid()]),input_id); - } - else _error_("could not recognize nature of vector from code " << code); - } - else if(M==iomodel->numberofelements+1){ - /*create transient input: */ - IssmDouble* times = xNew(N); - for(t=0;tSid()+t]; - switch(this->ObjectEnum()){ - case TriaEnum: transientinput->AddTimeInput(new TriaInput( input_enum,&value,P0Enum)); break; - case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,&value,P0Enum)); break; - case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,&value,P0Enum)); break; - default: _error_("Not implemented yet"); - } - } - datasetinput->AddInput(transientinput,input_id); - xDelete(times); - } - else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long"); - } - else if(vector_type==3){ //element vector - - /*For right now we are static */ - if(M==iomodel->numberofelements){ - /*create transient input: */ - IssmDouble* layers = xNewZeroInit(N);; - for(t=0;tSid()+t]; - DoubleArrayInput* arrayinput=new DoubleArrayInput(input_enum,layers,N); - datasetinput->AddInput(arrayinput,input_id); - xDelete(layers); - } - else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long"); - } - else _error_("Cannot add input for vector type " << vector_type << " (not supported)"); + case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P1Enum),input_id); break; + case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1Enum),input_id); break; + case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1Enum),input_id); break; + default: _error_("Not implemented yet"); + } + } + else if(N==this->GetNumberOfNodes(P0Enum) ){ + switch(this->ObjectEnum()){ + case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P0Enum),input_id); break; + case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P0Enum),input_id); break; + case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P0Enum),input_id); break; + default: _error_("Not implemented yet"); + } + } + else if(N==this->GetNumberOfNodes(P1xP2Enum)){ + switch(this->ObjectEnum()){ + case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P1xP2Enum),input_id); break; + case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1xP2Enum),input_id); break; + case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1xP2Enum),input_id); break; + default: _error_("Not implemented yet"); + } + } + else if(N==this->GetNumberOfNodes(P1xP3Enum)) { + switch(this->ObjectEnum()){ + case TriaEnum: datasetinput->AddInput(new TriaInput(input_enum,values,P1xP3Enum),input_id); break; + case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1xP3Enum),input_id); break; + case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1xP3Enum),input_id); break; + default: _error_("Not implemented yet"); + } + } + else _error_("Patch interpolation not supported yet"); + + } + else{ + _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long"); + } + + xDelete(values); + xDelete(vertexids); + } + else if(vector_type==2){ //element vector + + IssmDouble value; + + /*Are we in transient or static? */ + if(M==iomodel->numberofelements){ + if (code==5){ //boolean + datasetinput->AddInput(new BoolInput(input_enum,reCast(vector[this->Sid()])),input_id); + } + else if (code==6){ //integer + datasetinput->AddInput(new IntInput(input_enum,reCast(vector[this->Sid()])),input_id); + } + else if (code==7){ //IssmDouble + datasetinput->AddInput(new DoubleInput(input_enum,vector[this->Sid()]),input_id); + } + else _error_("could not recognize nature of vector from code " << code); + } + else if(M==iomodel->numberofelements+1){ + /*create transient input: */ + IssmDouble* times = xNew(N); + for(t=0;tSid()+t]; + switch(this->ObjectEnum()){ + case TriaEnum: transientinput->AddTimeInput(new TriaInput( input_enum,&value,P0Enum)); break; + case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,&value,P0Enum)); break; + case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,&value,P0Enum)); break; + default: _error_("Not implemented yet"); + } + } + datasetinput->AddInput(transientinput,input_id); + xDelete(times); + } + else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long"); + } + else if(vector_type==3){ //element vector + + /*For right now we are static */ + if(M==iomodel->numberofelements){ + /*create transient input: */ + IssmDouble* layers = xNewZeroInit(N);; + for(t=0;tSid()+t]; + DoubleArrayInput* arrayinput=new DoubleArrayInput(input_enum,layers,N); + datasetinput->AddInput(arrayinput,input_id); + xDelete(layers); + } + else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long"); + } + else _error_("Cannot add input for vector type " << vector_type << " (not supported)"); } /*}}}*/ @@ -3252,6 +3252,6 @@ /*Some intputs need to be computed, even if they are already in inputs, they might not be up to date!*/ switch(output_enum){ - case ViscousHeatingEnum: this->ViscousHeatingCreateInput(); break; - case StressMaxPrincipalEnum: this->StressMaxPrincipalCreateInput(); break; + case ViscousHeatingEnum: this->ViscousHeatingCreateInput(); break; + case StressMaxPrincipalEnum: this->StressMaxPrincipalCreateInput(); break; case StressTensorxxEnum: case StressTensorxyEnum: @@ -3266,5 +3266,5 @@ case StrainRateyzEnum: case StrainRatezzEnum: - case StrainRateeffectiveEnum: this->ComputeStrainRate(); break; + case StrainRateeffectiveEnum: this->ComputeStrainRate(); break; case DeviatoricStressxxEnum: case DeviatoricStressxyEnum: @@ -3287,26 +3287,26 @@ case CalvingrateyEnum: case CalvingCalvingrateEnum: - this->StrainRateparallel(); - this->StrainRateperpendicular(); - int calvinglaw; - this->FindParam(&calvinglaw,CalvingLawEnum); - switch(calvinglaw){ - case DefaultCalvingEnum: - //do nothing - break; - case CalvingLevermannEnum: - this->CalvingRateLevermann(); - break; - case CalvingVonmisesEnum: - case CalvingDev2Enum: - this->CalvingRateVonmises(); - break; - case CalvingCrevasseDepthEnum: - this->CalvingCrevasseDepth(); - break; - default: - _error_("Calving law "<StrainRateparallel(); + this->StrainRateperpendicular(); + int calvinglaw; + this->FindParam(&calvinglaw,CalvingLawEnum); + switch(calvinglaw){ + case DefaultCalvingEnum: + //do nothing + break; + case CalvingLevermannEnum: + this->CalvingRateLevermann(); + break; + case CalvingVonmisesEnum: + case CalvingDev2Enum: + this->CalvingRateVonmises(); + break; + case CalvingCrevasseDepthEnum: + this->CalvingCrevasseDepth(); + break; + default: + _error_("Calving law "<CalvingFluxLevelset(); break; case CalvingMeltingFluxLevelsetEnum: this->CalvingMeltingFluxLevelset(); break; @@ -3351,42 +3351,42 @@ switch(input->GetResultInterpolation()){ case P0Enum:{ - IssmDouble value; - bool bvalue; - Input* input = this->GetInput(output_enum); _assert_(input); - switch(input->ObjectEnum()){ - case DoubleInputEnum: - input->GetInputValue(&value); - break; - case BoolInputEnum: - input->GetInputValue(&bvalue); - value=reCast(bvalue); - break; - default: - Gauss* gauss = this->NewGauss(); - input->GetInputValue(&value,gauss); - delete gauss; - } - vector->SetValue(this->Sid(),value,INS_VAL); - break; - } + IssmDouble value; + bool bvalue; + Input* input = this->GetInput(output_enum); _assert_(input); + switch(input->ObjectEnum()){ + case DoubleInputEnum: + input->GetInputValue(&value); + break; + case BoolInputEnum: + input->GetInputValue(&bvalue); + value=reCast(bvalue); + break; + default: + Gauss* gauss = this->NewGauss(); + input->GetInputValue(&value,gauss); + delete gauss; + } + vector->SetValue(this->Sid(),value,INS_VAL); + break; + } case P1Enum:{ - const int NUM_VERTICES = this->GetNumberOfVertices(); - - IssmDouble *values = xNew(NUM_VERTICES); - int *connectivity= xNew(NUM_VERTICES); - int *sidlist = xNew(NUM_VERTICES); - - this->GetVerticesSidList(sidlist); - this->GetVerticesConnectivityList(connectivity); - this->GetInputListOnVertices(values,output_enum); - for(int i=0;i(connectivity[i]); - - vector->SetValues(NUM_VERTICES,sidlist,values,ADD_VAL); - - xDelete(values); - xDelete(connectivity); - xDelete(sidlist); - break; - } + const int NUM_VERTICES = this->GetNumberOfVertices(); + + IssmDouble *values = xNew(NUM_VERTICES); + int *connectivity= xNew(NUM_VERTICES); + int *sidlist = xNew(NUM_VERTICES); + + this->GetVerticesSidList(sidlist); + this->GetVerticesConnectivityList(connectivity); + this->GetInputListOnVertices(values,output_enum); + for(int i=0;i(connectivity[i]); + + vector->SetValues(NUM_VERTICES,sidlist,values,ADD_VAL); + + xDelete(values); + xDelete(connectivity); + xDelete(sidlist); + break; + } default: _error_("interpolation "<GetResultInterpolation())<<" not supported yet"); Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 24090) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 24091) @@ -243,8 +243,18 @@ /*Build B for this segment*/ - B[0] = +basis[index1]; - B[1] = +basis[index2]; - B[2] = -basis[index1]; - B[3] = -basis[index2]; + switch(finiteelement){ + case P0DGEnum: + B[0] = +basis[0]; + B[1] = -basis[0]; + break; + case P1DGEnum: + B[0] = +basis[index1]; + B[1] = +basis[index2]; + B[2] = -basis[index1]; + B[3] = -basis[index2]; + break; + default: + _error_("not supported yet"); + } /*Clean-up*/ @@ -268,8 +278,19 @@ /*Build B'*/ - Bprime[0] = basis[index1]; - Bprime[1] = basis[index2]; - Bprime[2] = basis[index1]; - Bprime[3] = basis[index2]; + /*Build B for this segment*/ + switch(finiteelement){ + case P0DGEnum: + Bprime[0] = basis[0]; + Bprime[1] = basis[0]; + break; + case P1DGEnum: + Bprime[0] = basis[index1]; + Bprime[1] = basis[index2]; + Bprime[2] = basis[index1]; + Bprime[3] = basis[index2]; + break; + default: + _error_("not supported yet"); + } /*Clean-up*/ @@ -305,4 +326,8 @@ switch(finiteelement){ + case P0Enum: case P0DGEnum: + basis[0]=triabasis[0]; + xDelete(triabasis); + return; case P1Enum: case P1DGEnum: basis[0]=triabasis[index1]; Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp (revision 24090) +++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp (revision 24091) @@ -96,5 +96,5 @@ int TriaInput::GetResultInterpolation(void){/*{{{*/ - if(this->interpolation_type==P0Enum){ + if(this->interpolation_type==P0Enum || this->interpolation_type==P0DGEnum){ return P0Enum; } Index: /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp (revision 24090) +++ /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp (revision 24091) @@ -62,5 +62,4 @@ /*FIXME: hardcode element degree for now*/ this->flux_degree= P1DGEnum; - //printf("-------------- file: Numericalflux.cpp line: %i\n",__LINE__); //this->flux_degree= P0DGEnum; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 24090) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 24091) @@ -512,35 +512,56 @@ } /*}}}*/ - if(finite_element==P1DGEnum && analysis!=UzawaPressureAnalysisEnum){/*Special case for DG...{{{*/ - int node_list[4]; - if(iomodel->domaintype!=Domain2DhorizontalEnum) _error_("not implemented yet"); - CreateEdges(iomodel); - CreateFaces(iomodel); - for(int i=0;inumberoffaces;i++){ - int e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2] - int e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2] - if(e2!=-2){ - if(epart[e1]!=epart[e2]){ - int i1=iomodel->faces[4*i+0]; - int i2=iomodel->faces[4*i+1]; - int pos=-1; - for(int j=0;j<3;j++) if(iomodel->elements[3*e2+j]==i1) pos=j; - if( pos==0){ node_list[0] = e2*3+0; node_list[1] = e2*3+2;} - else if(pos==1){ node_list[0] = e2*3+1; node_list[1] = e2*3+0;} - else if(pos==2){ node_list[0] = e2*3+2; node_list[1] = e2*3+1;} - else _error_("not supposed to happen"); - pos=-1; - for(int j=0;j<3;j++) if(iomodel->elements[3*e1+j]==i1) pos=j; - if( pos==0){ node_list[2] = e1*3+0; node_list[3] = e1*3+1;} - else if(pos==1){ node_list[2] = e1*3+1; node_list[3] = e1*3+2;} - else if(pos==2){ node_list[2] = e1*3+2; node_list[3] = e1*3+0;} - else _error_("not supposed to happen"); - for(int j=0;j<4;j++){ - int nid = node_list[j]; - AddNodeToRank(nodes_ranks,nodes_proc_count,nid,epart[e1]); - AddNodeToRank(nodes_ranks,nodes_proc_count,nid,epart[e2]); + if((finite_element==P0DGEnum || finite_element==P1DGEnum) && analysis!=UzawaPressureAnalysisEnum){/*Special case for DG...{{{*/ + if(finite_element==P1DGEnum){ + int node_list[4]; + if(iomodel->domaintype!=Domain2DhorizontalEnum) _error_("not implemented yet"); + CreateEdges(iomodel); + CreateFaces(iomodel); + for(int i=0;inumberoffaces;i++){ + int e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2] + int e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2] + if(e2!=-2){ + if(epart[e1]!=epart[e2]){ + int i1=iomodel->faces[4*i+0]; + int i2=iomodel->faces[4*i+1]; + int pos=-1; + for(int j=0;j<3;j++) if(iomodel->elements[3*e2+j]==i1) pos=j; + if( pos==0){ node_list[0] = e2*3+0; node_list[1] = e2*3+2;} + else if(pos==1){ node_list[0] = e2*3+1; node_list[1] = e2*3+0;} + else if(pos==2){ node_list[0] = e2*3+2; node_list[1] = e2*3+1;} + else _error_("not supposed to happen"); + pos=-1; + for(int j=0;j<3;j++) if(iomodel->elements[3*e1+j]==i1) pos=j; + if( pos==0){ node_list[2] = e1*3+0; node_list[3] = e1*3+1;} + else if(pos==1){ node_list[2] = e1*3+1; node_list[3] = e1*3+2;} + else if(pos==2){ node_list[2] = e1*3+2; node_list[3] = e1*3+0;} + else _error_("not supposed to happen"); + for(int j=0;j<4;j++){ + int nid = node_list[j]; + AddNodeToRank(nodes_ranks,nodes_proc_count,nid,epart[e1]); + AddNodeToRank(nodes_ranks,nodes_proc_count,nid,epart[e2]); + } } } } + } + else if(finite_element==P0DGEnum){ + int node_list[2]; + if(iomodel->domaintype!=Domain2DhorizontalEnum) _error_("not implemented yet"); + CreateEdges(iomodel); + CreateFaces(iomodel); + for(int i=0;inumberoffaces;i++){ + int e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2] + int e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2] + if(e2!=-2){ + if(epart[e1]!=epart[e2]){ + AddNodeToRank(nodes_ranks,nodes_proc_count,e2,epart[e1]); + //AddNodeToRank(nodes_ranks,nodes_proc_count,e1,epart[e2]); + } + } + } + } + else{ + _error_("not supported"); } }/*}}}*/