Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.cpp (revision 21310) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.cpp (revision 21311) @@ -963,4 +963,27 @@ } /*}}}*/ +void Element::GetInputListOnVerticesAtTime(IssmDouble* pvalue, int enumtype, IssmDouble time){/*{{{*/ + + /*Recover input*/ + Input* input=this->GetInput(enumtype); + if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); + + /*Fetch number vertices for this element*/ + int numvertices = this->GetNumberOfVertices(); + + /*Checks in debugging mode*/ + _assert_(pvalue); + + /* Start looping on the number of vertices: */ + Gauss*gauss=this->NewGauss(); + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(&pvalue[iv],gauss,time); + } + + /*clean-up*/ + delete gauss; +} +/*}}}*/ void Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/ @@ -1172,4 +1195,31 @@ this->GetInputListOnNodes(values,input_enum); vector->SetValues(numnodes,doflist,values,INS_VAL); + break; + default: + _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet"); + } + + /*Clean up*/ + xDelete(doflist); + xDelete(values); + +} +/*}}}*/ +void Element::GetVectorFromInputs(Vector* vector,int input_enum,int type, IssmDouble time){/*{{{*/ + + /*Fetch number vertices for this element and allocate arrays*/ + int numvertices = this->GetNumberOfVertices(); + int numnodes = this->GetNumberOfNodes(); + int* doflist = NULL; + IssmDouble* values = NULL; + + switch(type){ + case VertexPIdEnum: + doflist = xNew(numvertices); + values = xNew(numvertices); + /*Fill in values*/ + this->GetVertexPidList(doflist); + this->GetInputListOnVerticesAtTime(values,input_enum,time); + vector->SetValues(numvertices,doflist,values,INS_VAL); break; default: Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.h (revision 21310) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.h (revision 21311) @@ -89,4 +89,5 @@ void GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype); + void GetInputListOnVerticesAtTime(IssmDouble* pvalue,int enumtype,IssmDouble time); void GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue); void GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,IssmDouble* ug); @@ -100,4 +101,5 @@ void GetPhi(IssmDouble* phi, IssmDouble* epsilon, IssmDouble viscosity); void GetVectorFromInputs(Vector* vector, int name_enum, int type); + void GetVectorFromInputs(Vector* vector, int name_enum, int type,IssmDouble time); void GetVertexPidList(int* pidlist); void GetVerticesConnectivityList(int* connectivitylist); Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h (revision 21310) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h (revision 21311) @@ -64,4 +64,5 @@ void TransferForcing(FemModel* femmodel,int forcingenum); void TransferSealevel(FemModel* femmodel,int forcingenum); +void EarthMassTransport(FemModel* femmodel); //solution configuration Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp (revision 21310) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp (revision 21311) @@ -40,5 +40,5 @@ femmodel->parameters->FindParam(&isslr,TransientIsslrEnum); femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum); - + /*first, recover lat,long and radius vectors from vertices: */ VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); @@ -71,4 +71,14 @@ /*set configuration: */ if(isslr)femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum); + + /*figure out whether deltathickness was computed on the earth (mass trasnport module should + * have taken care of it for ice caps: */ + if(iscoupler){ + int ismasstransport,modelid,earthid; + femmodel->parameters->FindParam(&modelid,ModelIdEnum); + femmodel->parameters->FindParam(&earthid,EarthIdEnum); + femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum); + if((modelid==EarthIdEnum) & !ismasstransport) EarthMassTransport(femmodel); + } /*transfer deltathickness forcing from ice caps to earth model: */ @@ -359,2 +369,32 @@ } /*}}}*/ +void EarthMassTransport(FemModel* femmodel){ /*{{{*/ + + IssmDouble time,dt; + Vector *oldthickness = NULL; + Vector *newthickness = NULL; + Vector *deltathickness = NULL; + int nv; + + /*No mass transport module was called, so we are just going to retrieve the geometry thickness + * at this time step, at prior time step, and plug the difference as deltathickness: */ + femmodel->parameters->FindParam(&time,TimeEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + nv=femmodel->vertices->NumberOfVertices(); + + GetVectorFromInputsx(&newthickness,femmodel,ThicknessEnum,VertexSIdEnum); + GetVectorFromInputsx(&oldthickness,femmodel,ThicknessEnum,VertexSIdEnum,time-dt); + + /*compute deltathickness: */ + deltathickness = new Vector(nv); + newthickness->Copy(deltathickness); deltathickness->AXPY(oldthickness,-1); + + /*plug into elements:*/ + InputUpdateFromVectorx(femmodel,deltathickness,SealevelriseDeltathicknessEnum,VertexSIdEnum); + + /*free ressources:*/ + delete oldthickness; + delete newthickness; + delete deltathickness; + +} /*}}}*/ Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 21310) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp (revision 21311) @@ -33,4 +33,30 @@ *pvector=vector; } +void GetVectorFromInputsx(Vector** pvector,FemModel* femmodel,int name,int type,IssmDouble time){ + + int i; + Vector* vector=NULL; + + switch(type){ + case VertexPIdEnum: case VertexSIdEnum: + vector=new Vector(femmodel->vertices->NumberOfVertices()); + break; + case NodesEnum:case NodeSIdEnum: + vector=new Vector(femmodel->nodes->NumberOfNodes()); + break; + default: + _error_("vector type: " << EnumToStringx(type) << " not supported yet!"); + } + /*Look up in elements*/ + for(i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->GetVectorFromInputs(vector,name,type,time); + } + + vector->Assemble(); + + /*Assign output pointers:*/ + *pvector=vector; +} void GetVectorFromInputsx(IssmDouble** pvector,FemModel* femmodel,int name, int type){ Index: /issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h =================================================================== --- /issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h (revision 21310) +++ /issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h (revision 21311) @@ -9,4 +9,5 @@ /* local prototypes: */ void GetVectorFromInputsx( Vector** pvector,FemModel* femmodel,int name,int type); +void GetVectorFromInputsx( Vector** pvector,FemModel* femmodel,int name,int type,IssmDouble time); void GetVectorFromInputsx( IssmDouble** pvector,FemModel* femmodel,int name,int type);