Index: /issm/trunk/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 6215) +++ /issm/trunk/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 6216) @@ -58,8 +58,8 @@ if(iomodel->my_elements[i]){ element=(Element*)elements->GetObjectByOffset(counter); - element->Update(i,iomodel); //we need i to index into elements. + element->InputUpdateFromIoModel(i,iomodel); //we need i to index into elements. material=(Material*)materials->GetObjectByOffset(counter); - material->Update(i,iomodel); //we need i to index into elements. + material->InputUpdateFromIoModel(i,iomodel); //we need i to index into elements. counter++; } Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 6215) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 6216) @@ -56,5 +56,4 @@ virtual void DeleteResults(void)=0; virtual void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type)=0; - virtual void Update(int index, IoModel* iomodel)=0; virtual void UpdateGeometry(void)=0; virtual void InputToResult(int enum_type,int step,double time)=0; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6215) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6216) @@ -470,4 +470,174 @@ void Penta::InputUpdateFromVectorDakota(bool* vector, int name, int type){ ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Penta::InputUpdateFromIoModel(int index,IoModel* iomodel) {{{1*/ +void Penta::InputUpdateFromIoModel(int index,IoModel* iomodel){ + + /*Intermediaries*/ + IssmInt i,j; + int penta_vertex_ids[6]; + double nodeinputs[6]; + + /*Checks if debuging*/ + /*{{{2*/ + ISSMASSERT(iomodel->elements); + /*}}}*/ + + /*Recover vertices ids needed to initialize inputs*/ + for(i=0;i<6;i++){ + penta_vertex_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab + } + + //add as many inputs per element as requested: + if (iomodel->thickness) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,nodeinputs)); + } + if (iomodel->surface) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,nodeinputs)); + } + if (iomodel->bed) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(BedEnum,nodeinputs)); + } + if (iomodel->drag_coefficient) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(DragCoefficientEnum,nodeinputs)); + + if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index])); + if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index])); + this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type)); + + } + if (iomodel->melting_rate) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->melting_rate[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(MeltingRateEnum,nodeinputs)); + } + if (iomodel->accumulation_rate) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(AccumulationRateEnum,nodeinputs)); + } + if (iomodel->geothermalflux) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(GeothermalFluxEnum,nodeinputs)); + } + if (iomodel->pressure) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->pressure[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(PressureEnum,nodeinputs)); + } + if (iomodel->temperature) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->temperature[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,nodeinputs)); + } + if (iomodel->dhdt) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->dhdt[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(DhDtEnum,nodeinputs)); + } + /*vx,vy and vz: */ + if (iomodel->vx) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs)); + this->inputs->AddInput(new PentaVertexInput(VxOldEnum,nodeinputs)); + if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVxEnum,nodeinputs)); + } + if (iomodel->vy) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs)); + this->inputs->AddInput(new PentaVertexInput(VyOldEnum,nodeinputs)); + if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVyEnum,nodeinputs)); + } + if (iomodel->vz) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs)); + this->inputs->AddInput(new PentaVertexInput(VzOldEnum,nodeinputs)); + if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVzEnum,nodeinputs)); + } + if (iomodel->vx_obs) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(VxObsEnum,nodeinputs)); + } + if (iomodel->vy_obs) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(VyObsEnum,nodeinputs)); + } + if (iomodel->vz_obs) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(VzObsEnum,nodeinputs)); + } + if (iomodel->weights) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->weights[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(WeightsEnum,nodeinputs)); + } + if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index])); + if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index])); + if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index])); + if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index])); + + /*Control Inputs*/ + if (iomodel->control_analysis && iomodel->control_type){ + for(i=0;inum_control_type;i++){ + switch((int)iomodel->control_type[i]){ + case DhDtEnum: + if (iomodel->dhdt){ + for(j=0;j<6;j++)nodeinputs[j]=iomodel->dhdt[penta_vertex_ids[j]-1]/iomodel->yts; + this->inputs->AddInput(new ControlInput(DhDtEnum,PentaVertexInputEnum,nodeinputs)); + } + break; + case VxEnum: + if (iomodel->vx){ + for(j=0;j<6;j++)nodeinputs[j]=iomodel->vx[penta_vertex_ids[j]-1]/iomodel->yts; + this->inputs->AddInput(new ControlInput(VxEnum,PentaVertexInputEnum,nodeinputs)); + } + break; + case VyEnum: + if (iomodel->vy){ + for(j=0;j<6;j++)nodeinputs[j]=iomodel->vy[penta_vertex_ids[j]-1]/iomodel->yts; + this->inputs->AddInput(new ControlInput(VyEnum,PentaVertexInputEnum,nodeinputs)); + } + break; + case DragCoefficientEnum: + if (iomodel->drag_coefficient){ + for(j=0;j<6;j++)nodeinputs[j]=iomodel->drag_coefficient[penta_vertex_ids[j]-1]; + this->inputs->AddInput(new ControlInput(DragCoefficientEnum,PentaVertexInputEnum,nodeinputs)); + } + break; + case RheologyBbarEnum: + /*Matice will take care of it*/ break; + default: + ISSMERROR("Control %s not implemented yet",EnumToString((int)iomodel->control_type[i])); + } + } + } + + //Need to know the type of approximation for this element + if(iomodel->elements_type){ + if (*(iomodel->elements_type+index)==MacAyealApproximationEnum){ + this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealApproximationEnum)); + } + else if (*(iomodel->elements_type+index)==PattynApproximationEnum){ + this->inputs->AddInput(new IntInput(ApproximationEnum,PattynApproximationEnum)); + } + else if (*(iomodel->elements_type+index)==MacAyealPattynApproximationEnum){ + this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealPattynApproximationEnum)); + } + else if (*(iomodel->elements_type+index)==HutterApproximationEnum){ + this->inputs->AddInput(new IntInput(ApproximationEnum,HutterApproximationEnum)); + } + else if (*(iomodel->elements_type+index)==StokesApproximationEnum){ + this->inputs->AddInput(new IntInput(ApproximationEnum,StokesApproximationEnum)); + } + else if (*(iomodel->elements_type+index)==PattynStokesApproximationEnum){ + this->inputs->AddInput(new IntInput(ApproximationEnum,PattynStokesApproximationEnum)); + } + else if (*(iomodel->elements_type+index)==NoneApproximationEnum){ + this->inputs->AddInput(new IntInput(ApproximationEnum,NoneApproximationEnum)); + } + else{ + ISSMERROR("Approximation type %s not supported yet",EnumToString((int)*(iomodel->elements_type+index))); + } + } + } /*}}}*/ @@ -1740,5 +1910,5 @@ /*Fill with IoModel*/ - this->Update(index,iomodel); + this->InputUpdateFromIoModel(index,iomodel); /*Defaults if not provided in iomodel*/ @@ -1817,174 +1987,4 @@ } /*}}}*/ -/*FUNCTION Penta::Update(int index,IoModel* iomodel) {{{1*/ -void Penta::Update(int index,IoModel* iomodel){ - - /*Intermediaries*/ - IssmInt i,j; - int penta_vertex_ids[6]; - double nodeinputs[6]; - - /*Checks if debuging*/ - /*{{{2*/ - ISSMASSERT(iomodel->elements); - /*}}}*/ - - /*Recover vertices ids needed to initialize inputs*/ - for(i=0;i<6;i++){ - penta_vertex_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab - } - - //add as many inputs per element as requested: - if (iomodel->thickness) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,nodeinputs)); - } - if (iomodel->surface) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,nodeinputs)); - } - if (iomodel->bed) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(BedEnum,nodeinputs)); - } - if (iomodel->drag_coefficient) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(DragCoefficientEnum,nodeinputs)); - - if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index])); - if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index])); - this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type)); - - } - if (iomodel->melting_rate) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->melting_rate[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(MeltingRateEnum,nodeinputs)); - } - if (iomodel->accumulation_rate) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(AccumulationRateEnum,nodeinputs)); - } - if (iomodel->geothermalflux) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(GeothermalFluxEnum,nodeinputs)); - } - if (iomodel->pressure) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->pressure[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(PressureEnum,nodeinputs)); - } - if (iomodel->temperature) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->temperature[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,nodeinputs)); - } - if (iomodel->dhdt) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->dhdt[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(DhDtEnum,nodeinputs)); - } - /*vx,vy and vz: */ - if (iomodel->vx) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs)); - this->inputs->AddInput(new PentaVertexInput(VxOldEnum,nodeinputs)); - if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVxEnum,nodeinputs)); - } - if (iomodel->vy) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs)); - this->inputs->AddInput(new PentaVertexInput(VyOldEnum,nodeinputs)); - if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVyEnum,nodeinputs)); - } - if (iomodel->vz) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs)); - this->inputs->AddInput(new PentaVertexInput(VzOldEnum,nodeinputs)); - if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVzEnum,nodeinputs)); - } - if (iomodel->vx_obs) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(VxObsEnum,nodeinputs)); - } - if (iomodel->vy_obs) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(VyObsEnum,nodeinputs)); - } - if (iomodel->vz_obs) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(VzObsEnum,nodeinputs)); - } - if (iomodel->weights) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->weights[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(WeightsEnum,nodeinputs)); - } - if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index])); - if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index])); - if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index])); - if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index])); - - /*Control Inputs*/ - if (iomodel->control_analysis && iomodel->control_type){ - for(i=0;inum_control_type;i++){ - switch((int)iomodel->control_type[i]){ - case DhDtEnum: - if (iomodel->dhdt){ - for(j=0;j<6;j++)nodeinputs[j]=iomodel->dhdt[penta_vertex_ids[j]-1]/iomodel->yts; - this->inputs->AddInput(new ControlInput(DhDtEnum,PentaVertexInputEnum,nodeinputs)); - } - break; - case VxEnum: - if (iomodel->vx){ - for(j=0;j<6;j++)nodeinputs[j]=iomodel->vx[penta_vertex_ids[j]-1]/iomodel->yts; - this->inputs->AddInput(new ControlInput(VxEnum,PentaVertexInputEnum,nodeinputs)); - } - break; - case VyEnum: - if (iomodel->vy){ - for(j=0;j<6;j++)nodeinputs[j]=iomodel->vy[penta_vertex_ids[j]-1]/iomodel->yts; - this->inputs->AddInput(new ControlInput(VyEnum,PentaVertexInputEnum,nodeinputs)); - } - break; - case DragCoefficientEnum: - if (iomodel->drag_coefficient){ - for(j=0;j<6;j++)nodeinputs[j]=iomodel->drag_coefficient[penta_vertex_ids[j]-1]; - this->inputs->AddInput(new ControlInput(DragCoefficientEnum,PentaVertexInputEnum,nodeinputs)); - } - break; - case RheologyBbarEnum: - /*Matice will take care of it*/ break; - default: - ISSMERROR("Control %s not implemented yet",EnumToString((int)iomodel->control_type[i])); - } - } - } - - //Need to know the type of approximation for this element - if(iomodel->elements_type){ - if (*(iomodel->elements_type+index)==MacAyealApproximationEnum){ - this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealApproximationEnum)); - } - else if (*(iomodel->elements_type+index)==PattynApproximationEnum){ - this->inputs->AddInput(new IntInput(ApproximationEnum,PattynApproximationEnum)); - } - else if (*(iomodel->elements_type+index)==MacAyealPattynApproximationEnum){ - this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealPattynApproximationEnum)); - } - else if (*(iomodel->elements_type+index)==HutterApproximationEnum){ - this->inputs->AddInput(new IntInput(ApproximationEnum,HutterApproximationEnum)); - } - else if (*(iomodel->elements_type+index)==StokesApproximationEnum){ - this->inputs->AddInput(new IntInput(ApproximationEnum,StokesApproximationEnum)); - } - else if (*(iomodel->elements_type+index)==PattynStokesApproximationEnum){ - this->inputs->AddInput(new IntInput(ApproximationEnum,PattynStokesApproximationEnum)); - } - else if (*(iomodel->elements_type+index)==NoneApproximationEnum){ - this->inputs->AddInput(new IntInput(ApproximationEnum,NoneApproximationEnum)); - } - else{ - ISSMERROR("Approximation type %s not supported yet",EnumToString((int)*(iomodel->elements_type+index))); - } - } - -} -/*}}}*/ /*FUNCTION Penta::UpdateGeometry{{{1*/ void Penta::UpdateGeometry(void){ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 6215) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 6216) @@ -69,4 +69,5 @@ void InputUpdateFromVectorDakota(double* vector, int name, int type); void InputUpdateFromVectorDakota(int* vector, int name, int type); + void InputUpdateFromIoModel(int index, IoModel* iomodel); /*}}}*/ /*Element virtual functions definitions: {{{1*/ @@ -119,5 +120,4 @@ double SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); - void Update(int index, IoModel* iomodel); void UpdateGeometry(void); double TimeAdapt(); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6215) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6216) @@ -535,1702 +535,6 @@ } /*}}}*/ - -/*Element virtual functions definitions: */ -/*FUNCTION Tria::AverageOntoPartition {{{1*/ -void Tria::AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part){ - - bool already=false; - int i,j; - int partition[NUMVERTICES]; - int offset[NUMVERTICES]; - double area; - double mean; - double values[3]; - - /*First, get the area: */ - area=this->GetArea(); - - /*Figure out the average for this element: */ - this->GetDofList1(&offset[0]); - mean=0; - for(i=0;iFindParam(&analysis_counter,AnalysisCounterEnum); - - /*Get Element type*/ - this->element_type=this->element_type_list[analysis_counter]; - - /*Pick up nodes*/ - if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); - else this->nodes=NULL; - -} -/*}}}*/ -/*FUNCTION Tria::Configure {{{1*/ -void Tria::Configure(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ - - /*go into parameters and get the analysis_counter: */ - int analysis_counter; - parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); - - /*Get Element type*/ - this->element_type=this->element_type_list[analysis_counter]; - - /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective - * datasets, using internal ids and offsets hidden in hooks: */ - if(this->hnodes[analysis_counter]) this->hnodes[analysis_counter]->configure(nodesin); - this->hmatice->configure(materialsin); - this->hmatpar->configure(materialsin); - - /*Now, go pick up the objects inside the hooks: */ - if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); - else this->nodes=NULL; - this->matice=(Matice*)this->hmatice->delivers(); - this->matpar=(Matpar*)this->hmatpar->delivers(); - - /*point parameters to real dataset: */ - this->parameters=parametersin; - -} -/*}}}*/ -/*FUNCTION Tria::ControlInputGetGradient{{{1*/ -void Tria::ControlInputGetGradient(Vec gradient,int enum_type){ - - int doflist1[NUMVERTICES]; - Input* input=NULL; - - if(enum_type==RheologyBbarEnum){ - input=(Input*)matice->inputs->GetInput(enum_type); - } - else{ - input=inputs->GetInput(enum_type); - } - if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type)); - if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type)); - - this->GetDofList1(&doflist1[0]); - ((ControlInput*)input)->GetGradient(gradient,&doflist1[0]); - -}/*}}}*/ -/*FUNCTION Tria::ControlInputSetGradient{{{1*/ -void Tria::ControlInputSetGradient(double* gradient,int enum_type){ - - int doflist1[NUMVERTICES]; - double grad_list[NUMVERTICES]; - Input* grad_input=NULL; - Input* input=NULL; - - if(enum_type==RheologyBbarEnum){ - input=(Input*)matice->inputs->GetInput(enum_type); - } - else{ - input=inputs->GetInput(enum_type); - } - if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type)); - if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type)); - - this->GetDofList1(&doflist1[0]); - for(int i=0;iSetGradient(grad_input); - -}/*}}}*/ -/*FUNCTION Tria::RegularizeInversion {{{1*/ -double Tria::RegularizeInversion(){ - - /*constants: */ - const int numdof=NDOF2*NUMVERTICES; - - /* Intermediaries */ - int i,j,ig; - int num_controls; - int *control_type=NULL; - double Jelem = 0; - double cm_noisedmp; - double Jdet; - double xyz_list[NUMVERTICES][3]; - double dk[NDOF2],dB[NDOF2]; - GaussTria *gauss = NULL; - - /*retrieve parameters and inputs*/ - this->parameters->FindParam(&num_controls,NumControlsEnum); - this->parameters->FindParam(&control_type,NULL,ControlTypeEnum); - this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); - - /*If on water, return 0: */ - if(IsOnWater()) return 0; - - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - - for(i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - - /*Add Tikhonov regularization term to misfit: - * - * WARNING: for now, the regularization is only taken into account by the gradient - * the misfit reflects the response only! - * - * */ - if (control_type[i]==DragCoefficientEnum){ - Input* drag_input=inputs->GetInput(DragCoefficientEnum); ISSMASSERT(drag_input); - drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); - //Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight; - } - else if (control_type[i]==RheologyBbarEnum){ - //nothing - } - else if (control_type[i]==DhDtEnum){ - //nothing - } - else if (control_type[i]==VxEnum){ - //nothing - } - else if (control_type[i]==VyEnum){ - //nothing - } - else{ - ISSMERROR("unsupported control type: %s",EnumToString(control_type[i])); - } - } - - delete gauss; - } - - /*Clean up and return*/ - xfree((void**)&control_type); - return Jelem; -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrix {{{1*/ -void Tria::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){ - - /*retreive parameters: */ - ElementMatrix* Ke=NULL; - int analysis_type; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*Checks in debugging mode{{{2*/ - ISSMASSERT(this->nodes && this->matice && this->matpar && this->parameters && this->inputs); - /*}}}*/ - - /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - switch(analysis_type){ - case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum: - Ke=CreateKMatrixDiagnosticMacAyeal(); - break; - case DiagnosticHutterAnalysisEnum: - Ke=CreateKMatrixDiagnosticHutter(); - break; - case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: - Ke=CreateKMatrixSlope(); - break; - case PrognosticAnalysisEnum: - Ke=CreateKMatrixPrognostic(); - break; - case BalancedthicknessAnalysisEnum: - Ke=CreateKMatrixBalancedthickness(); - break; - case AdjointBalancedthicknessAnalysisEnum: - Ke=CreateKMatrixAdjointBalancedthickness(); - break; - case BalancedvelocitiesAnalysisEnum: - Ke=CreateKMatrixBalancedvelocities(); - break; - default: - ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); - } - - /*Add to global matrix*/ - if(Ke){ - Ke->AddToGlobal(Kgg,Kff,Kfs); - delete Ke; - } -} -/*}}}*/ -/*FUNCTION Tria::CreatePVector {{{1*/ -void Tria::CreatePVector(Vec pg, Vec pf){ - - /*retrive parameters: */ - ElementVector* pe=NULL; - int analysis_type; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*asserts: {{{*/ - /*if debugging mode, check that all pointers exist*/ - ISSMASSERT(this->nodes && this->matice && this->matpar && this->parameters && this->inputs); - /*}}}*/ - - /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ - switch(analysis_type){ - case DiagnosticHorizAnalysisEnum: - pe=CreatePVectorDiagnosticMacAyeal(); - break; - case AdjointHorizAnalysisEnum: - pe=CreatePVectorAdjointHoriz(); - break; - case DiagnosticHutterAnalysisEnum: - pe=CreatePVectorDiagnosticHutter(); - break; - case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: - pe=CreatePVectorSlope(); - break; - case PrognosticAnalysisEnum: - pe=CreatePVectorPrognostic(); - break; - case BalancedthicknessAnalysisEnum: - pe=CreatePVectorBalancedthickness(); - break; - case AdjointBalancedthicknessAnalysisEnum: - pe=CreatePVectorAdjointBalancedthickness(); - break; - case BalancedvelocitiesAnalysisEnum: - pe=CreatePVectorBalancedvelocities(); - break; - default: - ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); - } - - /*Add to global Vector*/ - if(pe){ - pe->AddToGlobal(pg,pf); - delete pe; - } -} -/*}}}*/ -/*FUNCTION Tria::DeleteResults {{{1*/ -void Tria::DeleteResults(void){ - - /*Delete and reinitialize results*/ - delete this->results; - this->results=new Results(); - -} -/*}}}*/ -/*FUNCTION Tria::GetNodeIndex {{{1*/ -int Tria::GetNodeIndex(Node* node){ - - ISSMASSERT(nodes); - for(int i=0;iGetParameterValue(&onbed,ElementOnBedEnum); - return onbed; -} -/*}}}*/ -/*FUNCTION Tria::IsOnShelf {{{1*/ -bool Tria::IsOnShelf(){ - - bool shelf; - inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); - return shelf; -} -/*}}}*/ -/*FUNCTION Tria::IsOnWater {{{1*/ -bool Tria::IsOnWater(){ - - bool water; - inputs->GetParameterValue(&water,ElementOnWaterEnum); - return water; -} -/*}}}*/ -/*FUNCTION Tria::GetSolutionFromInputs{{{1*/ -void Tria::GetSolutionFromInputs(Vec solution){ - - /*retrive parameters: */ - int analysis_type; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */ - if (analysis_type==DiagnosticHorizAnalysisEnum) - GetSolutionFromInputsDiagnosticHoriz(solution); - else if (analysis_type==DiagnosticHutterAnalysisEnum) - GetSolutionFromInputsDiagnosticHutter(solution); - else - ISSMERROR("analysis: %s not supported yet",EnumToString(analysis_type)); - -} -/*}}}*/ -/*FUNCTION Tria::GetVectorFromInputs{{{1*/ -void Tria::GetVectorFromInputs(Vec vector,int NameEnum){ - - int doflist1[NUMVERTICES]; - - /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */ - for(int i=0;iinputs->Size();i++){ - Input* input=(Input*)this->inputs->GetObjectByOffset(i); - if(input->EnumType()==NameEnum){ - /*We found the enum. Use its values to fill into the vector, using the vertices ids: */ - this->GetDofList1(&doflist1[0]); - input->GetVectorFromInputs(vector,&doflist1[0]); - break; - } - } -} -/*}}}*/ -/*FUNCTION Tria::Gradj {{{1*/ -void Tria::Gradj(Vec gradient,int control_type){ - - /*If on water, grad = 0: */ - if(IsOnWater())return; - - switch(control_type){ - case DragCoefficientEnum: - GradjDrag(gradient); - break; - case RheologyBbarEnum: - GradjB(gradient); - break; - case DhDtEnum: - GradjDhDt(gradient); - break; - case VxEnum: - GradjVx(gradient); - break; - case VyEnum: - GradjVy(gradient); - break; - default: - ISSMERROR("%s%i","control type not supported yet: ",control_type); - } -} -/*}}}*/ -/*FUNCTION Tria::GradjB{{{1*/ -void Tria::GradjB(Vec gradient){ - - /*Intermediaries*/ - int i,ig; - int doflist[NUMVERTICES]; - double vx,vy,lambda,mu,thickness,Jdet; - double cm_noisedmp,viscosity_complement; - double dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2]; - double xyz_list[NUMVERTICES][3]; - double basis[3],epsilon[3]; - double dbasis[NDOF2][NUMVERTICES]; - double grad_g[NUMVERTICES]; - double grad[NUMVERTICES]={0.0}; - GaussTria *gauss = NULL; - - /*retrieve some parameters: */ - this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList1(&doflist[0]); - - /*Retrieve all inputs*/ - Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); - Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); - Input* adjointx_input=inputs->GetInput(AdjointxEnum); ISSMASSERT(adjointx_input); - Input* adjointy_input=inputs->GetInput(AdjointyEnum); ISSMASSERT(adjointy_input); - Input* rheologyb_input=matice->inputs->GetInput(RheologyBbarEnum); ISSMASSERT(rheologyb_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(4); - for (ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - thickness_input->GetParameterValue(&thickness,gauss); - rheologyb_input->GetParameterDerivativeValue(&dB[0],&xyz_list[0][0],gauss); - vx_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); - vy_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); - adjointx_input->GetParameterDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss); - adjointy_input->GetParameterDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss); - - this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); - matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]); - - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(basis,gauss); - GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); - - /*standard gradient dJ/dki*/ - for (i=0;iweight*basis[i]; - } - /*Add regularization term*/ - for (i=0;iweight*(dbasis[0][i]*dB[0]+dbasis[1][i]*dB[1]); - for(i=0;iFindParam(&analysis_type,AnalysisTypeEnum); - - /*retrieve some parameters ands return if iceshelf: */ - this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); - if(IsOnShelf())return; - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList1(&doflist1[0]); - - /*Build frictoin element, needed later: */ - inputs->GetParameterValue(&drag_type,DragTypeEnum); - friction=new Friction("2d",inputs,matpar,analysis_type); - - /*Retrieve all inputs we will be needing: */ - Input* adjointx_input=inputs->GetInput(AdjointxEnum); ISSMASSERT(adjointx_input); - Input* adjointy_input=inputs->GetInput(AdjointyEnum); ISSMASSERT(adjointy_input); - Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); - Input* dragcoefficient_input=inputs->GetInput(DragCoefficientEnum); ISSMASSERT(dragcoefficient_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(4); - for (ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(l1l2l3, gauss); - GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss); - - /*Build alpha_complement_list: */ - if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); - else alpha_complement=0; - - dragcoefficient_input->GetParameterValue(&drag, gauss); - adjointx_input->GetParameterValue(&lambda, gauss); - adjointy_input->GetParameterValue(&mu, gauss); - vx_input->GetParameterValue(&vx,gauss); - vy_input->GetParameterValue(&vy,gauss); - dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); - - /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for (i=0;iweight*l1l2l3[i]; - - //noise dampening d/dki(1/2*(dk/dx)^2) - grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); - } - - /*Add gradje_g_gaussian vector to gradje_g: */ - for( i=0; iGetInput(AdjointEnum); ISSMASSERT(adjoint_input); - Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for (ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(l1l2l3, gauss); - - adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); - thickness_input->GetParameterValue(&thickness, gauss); - - for(i=0;iweight*l1l2l3[i]; - } - - VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); - - /*Clean up and return*/ - delete gauss; -} -/*}}}*/ -/*FUNCTION Tria::GradjVy{{{1*/ -void Tria::GradjVy(Vec gradient){ - - /*Intermediaries*/ - int i,ig; - int doflist1[NUMVERTICES]; - double thickness,Jdet; - double l1l2l3[3]; - double Dlambda[2]; - double xyz_list[NUMVERTICES][3]; - double grade_g[NUMVERTICES] = {0.0}; - GaussTria *gauss = NULL; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList1(&doflist1[0]); - - /*Retrieve all inputs we will be needing: */ - Input* adjoint_input=inputs->GetInput(AdjointEnum); ISSMASSERT(adjoint_input); - Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for (ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); - thickness_input->GetParameterValue(&thickness, gauss); - - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - GetNodalFunctions(l1l2l3, gauss); - - for(i=0;iweight*l1l2l3[i]; - } - - VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); - - /*Clean up and return*/ - delete gauss; -} -/*}}}*/ -/*FUNCTION Tria::InputControlUpdate{{{1*/ -void Tria::InputControlUpdate(double scalar,bool save_parameter){ - - /*Intermediary*/ - int num_controls; - int* control_type=NULL; - Input* input=NULL; - double cm_min,cm_max; - - /*retrieve some parameters: */ - this->parameters->FindParam(&cm_min,CmMinEnum); - this->parameters->FindParam(&cm_max,CmMaxEnum); - this->parameters->FindParam(&num_controls,NumControlsEnum); - this->parameters->FindParam(&control_type,NULL,ControlTypeEnum); - - for(int i=0;iinputs->GetInput(control_type[i]); ISSMASSERT(input); - } - else{ - input=(Input*)this->inputs->GetInput(control_type[i]); ISSMASSERT(input); - } - - if (input->Enum()!=ControlInputEnum){ - ISSMERROR("input %s is not a ControlInput",EnumToString(control_type[i])); - } - - ((ControlInput*)input)->UpdateValue(scalar); - input->Constrain(cm_min,cm_max); - if (save_parameter) ((ControlInput*)input)->SaveValue(); - - } - - /*Clean up and return*/ - xfree((void**)&control_type); -} -/*}}}*/ -/*FUNCTION Tria::InputConvergence{{{1*/ -bool Tria::InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){ - - bool converged=true; - int i; - Input** new_inputs=NULL; - Input** old_inputs=NULL; - - new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs - old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs - - for(i=0;iinputs->GetInput(enums[2*i+0]); - old_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+1]); - if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumToString(enums[2*i+0])); - if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumToString(enums[2*i+0])); - } - - /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/ - for(i=0;icriterionvalues[i]) converged=false; - } - - /*clean up and return*/ - xfree((void**)&new_inputs); - xfree((void**)&old_inputs); - return converged; -} -/*}}}*/ -/*FUNCTION Tria::InputDepthAverageAtBase {{{1*/ -void Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){ - - /*New input*/ - Input* oldinput=NULL; - Input* newinput=NULL; - - /*copy input of enum_type*/ - if (object_enum==ElementsEnum) - oldinput=(Input*)this->inputs->GetInput(enum_type); - else if (object_enum==MaterialsEnum) - oldinput=(Input*)this->matice->inputs->GetInput(enum_type); - else - ISSMERROR("object %s not supported yet",EnumToString(object_enum)); - if(!oldinput)ISSMERROR("%s%s"," could not find old input with enum: ",EnumToString(enum_type)); - newinput=(Input*)oldinput->copy(); - - /*Assign new name (average)*/ - newinput->ChangeEnum(average_enum_type); - - /*Add new input to current element*/ - if (object_enum==ElementsEnum) - this->inputs->AddInput((Input*)newinput); - else if (object_enum==MaterialsEnum) - this->matice->inputs->AddInput((Input*)newinput); - else - ISSMERROR("object %s not supported yet",EnumToString(object_enum)); -} -/*}}}*/ -/*FUNCTION Tria::InputDuplicate{{{1*/ -void Tria::InputDuplicate(int original_enum,int new_enum){ - - /*Call inputs method*/ - if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum); - -} -/*}}}*/ -/*FUNCTION Tria::InputScale{{{1*/ -void Tria::InputScale(int enum_type,double scale_factor){ - - Input* input=NULL; - - /*Make a copy of the original input: */ - input=(Input*)this->inputs->GetInput(enum_type); - if(!input)ISSMERROR(" could not find old input with enum: %s",EnumToString(enum_type)); - - /*Scale: */ - input->Scale(scale_factor); -} -/*}}}*/ -/*FUNCTION Tria::InputArtificialNoise{{{1*/ -void Tria::InputArtificialNoise(int enum_type,double min,double max){ - - Input* input=NULL; - - /*Make a copy of the original input: */ - input=(Input*)this->inputs->GetInput(enum_type); - if(!input)ISSMERROR(" could not find old input with enum: %s",EnumToString(enum_type)); - - /*ArtificialNoise: */ - input->ArtificialNoise(min,max); -} -/*}}}*/ -/*FUNCTION Tria::InputToResult{{{1*/ -void Tria::InputToResult(int enum_type,int step,double time){ - - int i; - Input *input = NULL; - - /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ - if (enum_type==RheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type); - else input=this->inputs->GetInput(enum_type); - if (!input) ISSMERROR("Input %s not found in tria->inputs",EnumToString(enum_type)); - - /*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result - * object out of the input, with the additional step and time information: */ - this->results->AddObject((Object*)input->SpawnResult(step,time)); - if(input->Enum()==ControlInputEnum) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time)); -} -/*}}}*/ -/*FUNCTION Tria::MassFlux {{{1*/ -double Tria::MassFlux( double* segment,bool process_units){ - - const int numdofs=2; - - int i; - double mass_flux=0; - double xyz_list[NUMVERTICES][3]; - double normal[2]; - double length,rho_ice; - double x1,y1,x2,y2,h1,h2; - double vx1,vx2,vy1,vy2; - GaussTria* gauss_1=NULL; - GaussTria* gauss_2=NULL; - - /*Get material parameters :*/ - rho_ice=matpar->GetRhoIce(); - - /*First off, check that this segment belongs to this element: */ - if ((int)*(segment+4)!=this->id)ISSMERROR("%s%i%s%i","error message: segment with id ",(int)*(segment+4)," does not belong to element with id:",this->id); - - /*Recover segment node locations: */ - x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3); - - /*Get xyz list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - - /*get area coordinates of 0 and 1 locations: */ - gauss_1=new GaussTria(); - gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]); - gauss_2=new GaussTria(); - gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]); - - normal[0]=cos(atan2(x1-x2,y2-y1)); - normal[1]=sin(atan2(x1-x2,y2-y1)); - - length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2)); - - Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); - Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); - - thickness_input->GetParameterValue(&h1, gauss_1); - thickness_input->GetParameterValue(&h2, gauss_2); - vx_input->GetParameterValue(&vx1,gauss_1); - vx_input->GetParameterValue(&vx2,gauss_2); - vy_input->GetParameterValue(&vy1,gauss_1); - vy_input->GetParameterValue(&vy2,gauss_2); - - mass_flux= rho_ice*length*( - (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+ - (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1] - ); - - /*Process units: */ - mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum,this->parameters); - - /*clean up and return:*/ - delete gauss_1; - delete gauss_2; - return mass_flux; -} -/*}}}*/ -/*FUNCTION Tria::MaxAbsVx{{{1*/ -void Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){ - - /*Get maximum:*/ - double maxabsvx=this->inputs->MaxAbs(VxEnum); - - /*process units if requested: */ - if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum,this->parameters); - - /*Assign output pointers:*/ - *pmaxabsvx=maxabsvx; -} -/*}}}*/ -/*FUNCTION Tria::MaxAbsVy{{{1*/ -void Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){ - - /*Get maximum:*/ - double maxabsvy=this->inputs->MaxAbs(VyEnum); - - /*process units if requested: */ - if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum,this->parameters); - - /*Assign output pointers:*/ - *pmaxabsvy=maxabsvy; -} -/*}}}*/ -/*FUNCTION Tria::MaxAbsVz{{{1*/ -void Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){ - - /*Get maximum:*/ - double maxabsvz=this->inputs->MaxAbs(VzEnum); - - /*process units if requested: */ - if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum,this->parameters); - - /*Assign output pointers:*/ - *pmaxabsvz=maxabsvz; -} -/*}}}*/ -/*FUNCTION Tria::MaxVel{{{1*/ -void Tria::MaxVel(double* pmaxvel, bool process_units){ - - /*Get maximum:*/ - double maxvel=this->inputs->Max(VelEnum); - - /*process units if requested: */ - if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum,this->parameters); - - /*Assign output pointers:*/ - *pmaxvel=maxvel; -} -/*}}}*/ -/*FUNCTION Tria::MaxVx{{{1*/ -void Tria::MaxVx(double* pmaxvx, bool process_units){ - - /*Get maximum:*/ - double maxvx=this->inputs->Max(VxEnum); - - /*process units if requested: */ - if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum,this->parameters); - - /*Assign output pointers:*/ - *pmaxvx=maxvx; -} -/*}}}*/ -/*FUNCTION Tria::MaxVy{{{1*/ -void Tria::MaxVy(double* pmaxvy, bool process_units){ - - /*Get maximum:*/ - double maxvy=this->inputs->Max(VyEnum); - - /*process units if requested: */ - if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum,this->parameters); - - /*Assign output pointers:*/ - *pmaxvy=maxvy; - -} -/*}}}*/ -/*FUNCTION Tria::MaxVz{{{1*/ -void Tria::MaxVz(double* pmaxvz, bool process_units){ - - /*Get maximum:*/ - double maxvz=this->inputs->Max(VzEnum); - - /*process units if requested: */ - if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum,this->parameters); - - /*Assign output pointers:*/ - *pmaxvz=maxvz; -} -/*}}}*/ -/*FUNCTION Tria::MinVel{{{1*/ -void Tria::MinVel(double* pminvel, bool process_units){ - - /*Get minimum:*/ - double minvel=this->inputs->Min(VelEnum); - - /*process units if requested: */ - if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum,this->parameters); - - /*Assign output pointers:*/ - *pminvel=minvel; -} -/*}}}*/ -/*FUNCTION Tria::MinVx{{{1*/ -void Tria::MinVx(double* pminvx, bool process_units){ - - /*Get minimum:*/ - double minvx=this->inputs->Min(VxEnum); - - /*process units if requested: */ - if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum,this->parameters); - - /*Assign output pointers:*/ - *pminvx=minvx; -} -/*}}}*/ -/*FUNCTION Tria::MinVy{{{1*/ -void Tria::MinVy(double* pminvy, bool process_units){ - - /*Get minimum:*/ - double minvy=this->inputs->Min(VyEnum); - - /*process units if requested: */ - if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum,this->parameters); - - /*Assign output pointers:*/ - *pminvy=minvy; -} -/*}}}*/ -/*FUNCTION Tria::MinVz{{{1*/ -void Tria::MinVz(double* pminvz, bool process_units){ - - /*Get minimum:*/ - double minvz=this->inputs->Min(VzEnum); - - /*process units if requested: */ - if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum,this->parameters); - - /*Assign output pointers:*/ - *pminvz=minvz; -} -/*}}}*/ -/*FUNCTION Tria::TimeAdapt{{{1*/ -double Tria::TimeAdapt(void){ - - /*intermediary: */ - int i; - double C,dt; - double dx,dy; - double maxx,minx; - double maxy,miny; - double maxabsvx,maxabsvy; - double xyz_list[NUMVERTICES][3]; - - /*get CFL coefficient:*/ - this->parameters->FindParam(&C,CflCoefficientEnum); - - /*Get for Vx and Vy, the max of abs value: */ - this->MaxAbsVx(&maxabsvx,false); - this->MaxAbsVy(&maxabsvy,false); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], this->nodes, NUMVERTICES); - - minx=xyz_list[0][0]; - maxx=xyz_list[0][0]; - miny=xyz_list[0][1]; - maxy=xyz_list[0][1]; - - for(i=1;imaxx)maxx=xyz_list[i][0]; - if (xyz_list[i][1]maxy)maxy=xyz_list[i][1]; - } - dx=maxx-minx; - dy=maxy-miny; - - /*CFL criterion: */ - dt=C/(maxabsvy/dx+maxabsvy/dy); - - return dt; -} -/*}}}*/ -/*FUNCTION Tria::ThicknessAbsMisfit {{{1*/ -double Tria::ThicknessAbsMisfit(bool process_units){ - - /* Constants */ - const int numdof=1*NUMVERTICES; - - /*Intermediaries*/ - int i,ig; - double thickness,thicknessobs,weight; - double Jdet; - double Jelem = 0; - double xyz_list[NUMVERTICES][3]; - GaussTria *gauss = NULL; - - /*If on water, return 0: */ - if(IsOnWater())return 0; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - - /*Retrieve all inputs we will be needing: */ - Input* thickness_input =inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); - Input* thicknessobs_input=inputs->GetInput(ThicknessObsEnum);ISSMASSERT(thicknessobs_input); - Input* weights_input =inputs->GetInput(WeightsEnum); ISSMASSERT(weights_input); - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for (ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - - /*Get parameters at gauss point*/ - thickness_input->GetParameterValue(&thickness,gauss); - thicknessobs_input->GetParameterValue(&thicknessobs,gauss); - weights_input->GetParameterValue(&weight,gauss); - - /*compute ThicknessAbsMisfit*/ - Jelem+=0.5*pow(thickness-thicknessobs,2.0)*weight*Jdet*gauss->weight; - } - - /* clean up and Return: */ - delete gauss; - return Jelem; -} -/*}}}*/ -/*FUNCTION Tria::SurfaceAbsVelMisfit {{{1*/ -double Tria::SurfaceAbsVelMisfit(bool process_units){ - - const int numdof=NDOF2*NUMVERTICES; - - int i,ig; - double Jelem=0; - double velocity_mag,obs_velocity_mag; - double meanvel, epsvel,misfit,Jdet; - double xyz_list[NUMVERTICES][3]; - double vx_list[NUMVERTICES]; - double vy_list[NUMVERTICES]; - double obs_vx_list[NUMVERTICES]; - double obs_vy_list[NUMVERTICES]; - double misfit_square_list[NUMVERTICES]; - double misfit_list[NUMVERTICES]; - double weights_list[NUMVERTICES]; - GaussTria *gauss=NULL; - - /*If on water, return 0: */ - if(IsOnWater())return 0; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - - /* Recover input data: */ - GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); - GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); - GetParameterListOnVertices(&vx_list[0],VxEnum); - GetParameterListOnVertices(&vy_list[0],VyEnum); - GetParameterListOnVertices(&weights_list[0],WeightsEnum); - - /*retrieve some parameters: */ - this->parameters->FindParam(&meanvel,MeanVelEnum); - this->parameters->FindParam(&epsvel,EpsVelEnum); - - /* Compute SurfaceAbsVelMisfit at the 3 nodes - * Here we integrate linearized functions: - * - * J(E) = int_E sum_{i=1}^3 J_i Phi_i - * - * where J_i are the misfits at the 3 nodes of the triangle - * Phi_i is the nodal function (P1) with respect to - * the vertex i - */ - - /*We are using an absolute misfit: - * - * 1 [ 2 2 ] - * J = --- | (u - u ) + (v - v ) | - * 2 [ obs obs ] - * - */ - for (i=0;iparameters); - - /*Apply weights to misfits*/ - for (i=0;ibegin();igend(); ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); - Jelem+=misfit*Jdet*gauss->weight; - } - - /*clean up and Return: */ - delete gauss; - return Jelem; -} -/*}}}*/ -/*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/ -double Tria::SurfaceRelVelMisfit(bool process_units){ - - const int numdof=2*NUMVERTICES; - - int i,ig; - double Jelem=0; - double scalex=1,scaley=1; - double meanvel, epsvel,misfit,Jdet; - double velocity_mag,obs_velocity_mag; - double xyz_list[NUMVERTICES][3]; - double vx_list[NUMVERTICES]; - double vy_list[NUMVERTICES]; - double obs_vx_list[NUMVERTICES]; - double obs_vy_list[NUMVERTICES]; - double misfit_square_list[NUMVERTICES]; - double misfit_list[NUMVERTICES]; - double weights_list[NUMVERTICES]; - GaussTria *gauss=NULL; - - /*If on water, return 0: */ - if(IsOnWater())return 0; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - - /* Recover input data: */ - GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); - GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); - GetParameterListOnVertices(&vx_list[0],VxEnum); - GetParameterListOnVertices(&vy_list[0],VyEnum); - GetParameterListOnVertices(&weights_list[0],WeightsEnum); - - /*retrieve some parameters: */ - this->parameters->FindParam(&meanvel,MeanVelEnum); - this->parameters->FindParam(&epsvel,EpsVelEnum); - - /* Compute SurfaceRelVelMisfit at the 3 nodes - * Here we integrate linearized functions: - * - * J(E) = int_E sum_{i=1}^3 J_i Phi_i - * - * where J_i are the misfits at the 3 nodes of the triangle - * Phi_i is the nodal function (P1) with respect to - * the vertex i - */ - - /*We are using a relative misfit: - * - * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] - * J = --- | ------------- (u - u ) + ------------- (v - v ) | - * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] - * obs obs - */ - for (i=0;iparameters); - - /*Apply weights to misfits*/ - for (i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); - Jelem+=misfit*Jdet*gauss->weight; - } - - /*clean up and Return: */ - delete gauss; - return Jelem; -} -/*}}}*/ -/*FUNCTION Tria::SurfaceLogVelMisfit {{{1*/ -double Tria::SurfaceLogVelMisfit(bool process_units){ - - const int numdof=NDOF2*NUMVERTICES; - - int i,ig; - double Jelem=0; - double scalex=1,scaley=1; - double meanvel, epsvel,misfit,Jdet; - double velocity_mag,obs_velocity_mag; - double xyz_list[NUMVERTICES][3]; - double vx_list[NUMVERTICES]; - double vy_list[NUMVERTICES]; - double obs_vx_list[NUMVERTICES]; - double obs_vy_list[NUMVERTICES]; - double misfit_square_list[NUMVERTICES]; - double misfit_list[NUMVERTICES]; - double weights_list[NUMVERTICES]; - GaussTria *gauss=NULL; - - /*If on water, return 0: */ - if(IsOnWater())return 0; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - - /* Recover input data: */ - GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); - GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); - GetParameterListOnVertices(&vx_list[0],VxEnum); - GetParameterListOnVertices(&vy_list[0],VyEnum); - GetParameterListOnVertices(&weights_list[0],WeightsEnum); - - /*retrieve some parameters: */ - this->parameters->FindParam(&meanvel,MeanVelEnum); - this->parameters->FindParam(&epsvel,EpsVelEnum); - - /* Compute SurfaceLogVelMisfit at the 3 nodes - * Here we integrate linearized functions: - * - * J(E) = int_E sum_{i=1}^3 J_i Phi_i - * - * where J_i are the misfits at the 3 nodes of the triangle - * Phi_i is the nodal function (P1) with respect to - * the vertex i - */ - - /*We are using a logarithmic misfit: - * - * [ vel + eps ] 2 - * J = 4 \bar{v}^2 | log ( ----------- ) | - * [ vel + eps ] - * obs - */ - for (i=0;iparameters); - - /*Apply weights to misfits*/ - for (i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); - Jelem+=misfit*Jdet*gauss->weight; - } - - /*clean-up and Return: */ - delete gauss; - return Jelem; -} -/*}}}*/ -/*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/ -double Tria::SurfaceLogVxVyMisfit(bool process_units){ - - const int numdof=NDOF2*NUMVERTICES; - - int i,ig; - int fit=-1; - double Jelem=0, S=0; - double scalex=1,scaley=1; - double meanvel, epsvel, misfit, Jdet; - double velocity_mag,obs_velocity_mag; - double xyz_list[NUMVERTICES][3]; - double vx_list[NUMVERTICES]; - double vy_list[NUMVERTICES]; - double obs_vx_list[NUMVERTICES]; - double obs_vy_list[NUMVERTICES]; - double misfit_square_list[NUMVERTICES]; - double misfit_list[NUMVERTICES]; - double weights_list[NUMVERTICES]; - GaussTria *gauss=NULL; - - /*If on water, return 0: */ - if(IsOnWater())return 0; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - - /* Recover input data: */ - GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); - GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); - GetParameterListOnVertices(&vx_list[0],VxEnum); - GetParameterListOnVertices(&vy_list[0],VyEnum); - GetParameterListOnVertices(&weights_list[0],WeightsEnum); - - /*retrieve some parameters: */ - this->parameters->FindParam(&meanvel,MeanVelEnum); - this->parameters->FindParam(&epsvel,EpsVelEnum); - - /* Compute SurfaceLogVxVyMisfit at the 3 nodes - * Here we integrate linearized functions: - * - * J(E) = int_E sum_{i=1}^3 J_i Phi_i - * - * where J_i are the misfits at the 3 nodes of the triangle - * Phi_i is the nodal function (P1) with respect to - * the vertex i - */ - - /*We are using an logarithmic 2 misfit: - * - * 1 [ |u| + eps 2 |v| + eps 2 ] - * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | - * 2 [ |u |+ eps |v |+ eps ] - * obs obs - */ - for (i=0;iparameters); - - /*Apply weights to misfits*/ - for (i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); - Jelem+=misfit*Jdet*gauss->weight; - } - - /*clean-up and Return: */ - delete gauss; - return Jelem; -} -/*}}}*/ -/*FUNCTION Tria::SurfaceAverageVelMisfit {{{1*/ -double Tria::SurfaceAverageVelMisfit(bool process_units){ - - const int numdof=2*NUMVERTICES; - - int i,ig; - int fit=-1; - double Jelem=0,S=0; - double scalex=1, scaley=1; - double meanvel, epsvel,Jdet; - double velocity_mag,obs_velocity_mag,misfit; - double xyz_list[NUMVERTICES][3]; - double vx_list[NUMVERTICES]; - double vy_list[NUMVERTICES]; - double obs_vx_list[NUMVERTICES]; - double obs_vy_list[NUMVERTICES]; - double misfit_square_list[NUMVERTICES]; - double misfit_list[NUMVERTICES]; - double weights_list[NUMVERTICES]; - GaussTria *gauss=NULL; - - /*If on water, return 0: */ - if(IsOnWater())return 0; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - - /* Recover input data: */ - inputs->GetParameterValue(&S,SurfaceAreaEnum); - GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); - GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); - GetParameterListOnVertices(&vx_list[0],VxEnum); - GetParameterListOnVertices(&vy_list[0],VyEnum); - GetParameterListOnVertices(&weights_list[0],WeightsEnum); - - /*retrieve some parameters: */ - this->parameters->FindParam(&meanvel,MeanVelEnum); - this->parameters->FindParam(&epsvel,EpsVelEnum); - - /* Compute SurfaceAverageVelMisfit at the 3 nodes - * Here we integrate linearized functions: - * - * J(E) = int_E sum_{i=1}^3 J_i Phi_i - * - * where J_i are the misfits at the 3 nodes of the triangle - * Phi_i is the nodal function (P1) with respect to - * the vertex i - */ - - /*We are using a spacially average absolute misfit: - * - * 1 2 2 - * J = --- sqrt( (u - u ) + (v - v ) ) - * S obs obs - */ - for (i=0;iparameters); - - /*Take the square root, and scale by surface: */ - for (i=0;ibegin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); - Jelem+=misfit*Jdet*gauss->weight; - } - - /*clean-up and Return: */ - delete gauss; - return Jelem; -} -/*}}}*/ -/*FUNCTION Tria::PatchFill{{{1*/ -void Tria::PatchFill(int* prow, Patch* patch){ - - int i,row; - int vertices_ids[3]; - - /*recover pointer: */ - row=*prow; - - for(i=0;i<3;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch. - - for(i=0;iresults->Size();i++){ - ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); - - /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand - *it to the result object, to fill the rest: */ - patch->fillelementinfo(row,this->id,vertices_ids,3); - elementresult->PatchFill(row,patch); - - /*increment rower: */ - row++; - } - - /*Assign output pointers:*/ - *prow=row; -} -/*}}}*/ -/*FUNCTION Tria::PatchSize{{{1*/ -void Tria::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){ - - int i; - int numrows = 0; - int numnodes = 0; - - /*Go through all the results objects, and update the counters: */ - for (i=0;iresults->Size();i++){ - ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); - /*first, we have one more result: */ - numrows++; - /*now, how many vertices and how many nodal values for this result? :*/ - numnodes=elementresult->NumberOfNodalValues(); //ask result object. - } - - /*Assign output pointers:*/ - *pnumrows=numrows; - *pnumvertices=NUMVERTICES; - *pnumnodes=numnodes; -} -/*}}}*/ -/*FUNCTION Tria::ProcessResultsUnits{{{1*/ -void Tria::ProcessResultsUnits(void){ - - int i; - - for(i=0;iresults->Size();i++){ - ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); - elementresult->ProcessUnits(this->parameters); - } -} -/*}}}*/ -/*FUNCTION Tria::SurfaceArea {{{1*/ -double Tria::SurfaceArea(void){ - - int i; - double S; - double normal[3]; - double v13[3],v23[3]; - double xyz_list[NUMVERTICES][3]; - - /*If on water, return 0: */ - if(IsOnWater())return 0; - - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - - for (i=0;i<3;i++){ - v13[i]=xyz_list[0][i]-xyz_list[2][i]; - v23[i]=xyz_list[1][i]-xyz_list[2][i]; - } - - normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; - normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; - normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; - - S = 0.5 * sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2)); - - /*Return: */ - return S; -} -/*}}}*/ -/*FUNCTION Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){{{1*/ -void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index - - /*Intermediaries*/ - int i,j; - int tria_node_ids[3]; - int tria_vertex_ids[3]; - int tria_type; - double nodeinputs[3]; - - /*Checks if debuging*/ - /*{{{2*/ - ISSMASSERT(iomodel->elements); - /*}}}*/ - - /*Recover element type*/ - if ((analysis_type==PrognosticAnalysisEnum || analysis_type==BalancedthicknessAnalysisEnum) && iomodel->prognostic_DG){ - /*P1 Discontinuous Galerkin*/ - tria_type=P1DGEnum; - } - else{ - /*P1 Continuous Galerkin*/ - tria_type=P1Enum; - } - this->SetElementType(tria_type,analysis_counter); - - /*Recover vertices ids needed to initialize inputs*/ - for(i=0;i<3;i++){ - tria_vertex_ids[i]=(int)iomodel->elements[3*index+i]; //ids for vertices are in the elements array from Matlab - } - - /*Recover nodes ids needed to initialize the node hook.*/ - if (tria_type==P1DGEnum){ - /*Discontinuous Galerkin*/ - tria_node_ids[0]=iomodel->nodecounter+3*index+1; - tria_node_ids[1]=iomodel->nodecounter+3*index+2; - tria_node_ids[2]=iomodel->nodecounter+3*index+3; - } - else{ - /*Continuous Galerkin*/ - for(i=0;i<3;i++){ - tria_node_ids[i]=iomodel->nodecounter+(int)*(iomodel->elements+3*index+i); //ids for vertices are in the elements array from Matlab - } - } - - /*hooks: */ - this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type - - /*Fill with IoModel*/ - this->Update(index,iomodel); - - /*Defaults if not provided in iomodel*/ - switch(analysis_type){ - - case DiagnosticHorizAnalysisEnum: - - /*default vx,vy and vz: either observation or 0 */ - if(!iomodel->vx){ - if (iomodel->vx_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts; - else for(i=0;i<3;i++)nodeinputs[i]=0; - this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs)); - this->inputs->AddInput(new TriaVertexInput(VxOldEnum,nodeinputs)); - if(iomodel->qmu_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVxEnum,nodeinputs)); - } - if(!iomodel->vy){ - if (iomodel->vy_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts; - else for(i=0;i<3;i++)nodeinputs[i]=0; - this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs)); - this->inputs->AddInput(new TriaVertexInput(VyOldEnum,nodeinputs)); - if(iomodel->qmu_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVyEnum,nodeinputs)); - } - if(!iomodel->vz){ - if (iomodel->vz_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_vertex_ids[i]-1]/iomodel->yts; - else for(i=0;i<3;i++)nodeinputs[i]=0; - this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs)); - this->inputs->AddInput(new TriaVertexInput(VzOldEnum,nodeinputs)); - if(iomodel->qmu_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVzEnum,nodeinputs)); - } - if(!iomodel->pressure){ - for(i=0;i<3;i++)nodeinputs[i]=0; - if(iomodel->qmu_analysis){ - this->inputs->AddInput(new TriaVertexInput(PressureEnum,nodeinputs)); - this->inputs->AddInput(new TriaVertexInput(QmuPressureEnum,nodeinputs)); - } - } - break; - - default: - /*No update for other solution types*/ - break; - - } - - //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. - this->parameters=NULL; -} -/*}}}*/ -/*FUNCTION Tria::Update(int index, IoModel* iomodel){{{1*/ -void Tria::Update(int index, IoModel* iomodel){ //i is the element index +/*FUNCTION Tria::InputUpdateFromIoModel(int index, IoModel* iomodel){{{1*/ +void Tria::InputUpdateFromIoModel(int index, IoModel* iomodel){ //i is the element index /*Intermediaries*/ @@ -2375,4 +679,1700 @@ } /*}}}*/ + +/*Element virtual functions definitions: */ +/*FUNCTION Tria::AverageOntoPartition {{{1*/ +void Tria::AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part){ + + bool already=false; + int i,j; + int partition[NUMVERTICES]; + int offset[NUMVERTICES]; + double area; + double mean; + double values[3]; + + /*First, get the area: */ + area=this->GetArea(); + + /*Figure out the average for this element: */ + this->GetDofList1(&offset[0]); + mean=0; + for(i=0;iFindParam(&analysis_counter,AnalysisCounterEnum); + + /*Get Element type*/ + this->element_type=this->element_type_list[analysis_counter]; + + /*Pick up nodes*/ + if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); + else this->nodes=NULL; + +} +/*}}}*/ +/*FUNCTION Tria::Configure {{{1*/ +void Tria::Configure(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ + + /*go into parameters and get the analysis_counter: */ + int analysis_counter; + parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); + + /*Get Element type*/ + this->element_type=this->element_type_list[analysis_counter]; + + /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective + * datasets, using internal ids and offsets hidden in hooks: */ + if(this->hnodes[analysis_counter]) this->hnodes[analysis_counter]->configure(nodesin); + this->hmatice->configure(materialsin); + this->hmatpar->configure(materialsin); + + /*Now, go pick up the objects inside the hooks: */ + if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); + else this->nodes=NULL; + this->matice=(Matice*)this->hmatice->delivers(); + this->matpar=(Matpar*)this->hmatpar->delivers(); + + /*point parameters to real dataset: */ + this->parameters=parametersin; + +} +/*}}}*/ +/*FUNCTION Tria::ControlInputGetGradient{{{1*/ +void Tria::ControlInputGetGradient(Vec gradient,int enum_type){ + + int doflist1[NUMVERTICES]; + Input* input=NULL; + + if(enum_type==RheologyBbarEnum){ + input=(Input*)matice->inputs->GetInput(enum_type); + } + else{ + input=inputs->GetInput(enum_type); + } + if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type)); + if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type)); + + this->GetDofList1(&doflist1[0]); + ((ControlInput*)input)->GetGradient(gradient,&doflist1[0]); + +}/*}}}*/ +/*FUNCTION Tria::ControlInputSetGradient{{{1*/ +void Tria::ControlInputSetGradient(double* gradient,int enum_type){ + + int doflist1[NUMVERTICES]; + double grad_list[NUMVERTICES]; + Input* grad_input=NULL; + Input* input=NULL; + + if(enum_type==RheologyBbarEnum){ + input=(Input*)matice->inputs->GetInput(enum_type); + } + else{ + input=inputs->GetInput(enum_type); + } + if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type)); + if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type)); + + this->GetDofList1(&doflist1[0]); + for(int i=0;iSetGradient(grad_input); + +}/*}}}*/ +/*FUNCTION Tria::RegularizeInversion {{{1*/ +double Tria::RegularizeInversion(){ + + /*constants: */ + const int numdof=NDOF2*NUMVERTICES; + + /* Intermediaries */ + int i,j,ig; + int num_controls; + int *control_type=NULL; + double Jelem = 0; + double cm_noisedmp; + double Jdet; + double xyz_list[NUMVERTICES][3]; + double dk[NDOF2],dB[NDOF2]; + GaussTria *gauss = NULL; + + /*retrieve parameters and inputs*/ + this->parameters->FindParam(&num_controls,NumControlsEnum); + this->parameters->FindParam(&control_type,NULL,ControlTypeEnum); + this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); + + /*If on water, return 0: */ + if(IsOnWater()) return 0; + + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + + for(i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + + /*Add Tikhonov regularization term to misfit: + * + * WARNING: for now, the regularization is only taken into account by the gradient + * the misfit reflects the response only! + * + * */ + if (control_type[i]==DragCoefficientEnum){ + Input* drag_input=inputs->GetInput(DragCoefficientEnum); ISSMASSERT(drag_input); + drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); + //Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight; + } + else if (control_type[i]==RheologyBbarEnum){ + //nothing + } + else if (control_type[i]==DhDtEnum){ + //nothing + } + else if (control_type[i]==VxEnum){ + //nothing + } + else if (control_type[i]==VyEnum){ + //nothing + } + else{ + ISSMERROR("unsupported control type: %s",EnumToString(control_type[i])); + } + } + + delete gauss; + } + + /*Clean up and return*/ + xfree((void**)&control_type); + return Jelem; +} +/*}}}*/ +/*FUNCTION Tria::CreateKMatrix {{{1*/ +void Tria::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){ + + /*retreive parameters: */ + ElementMatrix* Ke=NULL; + int analysis_type; + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*Checks in debugging mode{{{2*/ + ISSMASSERT(this->nodes && this->matice && this->matpar && this->parameters && this->inputs); + /*}}}*/ + + /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ + switch(analysis_type){ + case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum: + Ke=CreateKMatrixDiagnosticMacAyeal(); + break; + case DiagnosticHutterAnalysisEnum: + Ke=CreateKMatrixDiagnosticHutter(); + break; + case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: + Ke=CreateKMatrixSlope(); + break; + case PrognosticAnalysisEnum: + Ke=CreateKMatrixPrognostic(); + break; + case BalancedthicknessAnalysisEnum: + Ke=CreateKMatrixBalancedthickness(); + break; + case AdjointBalancedthicknessAnalysisEnum: + Ke=CreateKMatrixAdjointBalancedthickness(); + break; + case BalancedvelocitiesAnalysisEnum: + Ke=CreateKMatrixBalancedvelocities(); + break; + default: + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); + } + + /*Add to global matrix*/ + if(Ke){ + Ke->AddToGlobal(Kgg,Kff,Kfs); + delete Ke; + } +} +/*}}}*/ +/*FUNCTION Tria::CreatePVector {{{1*/ +void Tria::CreatePVector(Vec pg, Vec pf){ + + /*retrive parameters: */ + ElementVector* pe=NULL; + int analysis_type; + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*asserts: {{{*/ + /*if debugging mode, check that all pointers exist*/ + ISSMASSERT(this->nodes && this->matice && this->matpar && this->parameters && this->inputs); + /*}}}*/ + + /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ + switch(analysis_type){ + case DiagnosticHorizAnalysisEnum: + pe=CreatePVectorDiagnosticMacAyeal(); + break; + case AdjointHorizAnalysisEnum: + pe=CreatePVectorAdjointHoriz(); + break; + case DiagnosticHutterAnalysisEnum: + pe=CreatePVectorDiagnosticHutter(); + break; + case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum: + pe=CreatePVectorSlope(); + break; + case PrognosticAnalysisEnum: + pe=CreatePVectorPrognostic(); + break; + case BalancedthicknessAnalysisEnum: + pe=CreatePVectorBalancedthickness(); + break; + case AdjointBalancedthicknessAnalysisEnum: + pe=CreatePVectorAdjointBalancedthickness(); + break; + case BalancedvelocitiesAnalysisEnum: + pe=CreatePVectorBalancedvelocities(); + break; + default: + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type)); + } + + /*Add to global Vector*/ + if(pe){ + pe->AddToGlobal(pg,pf); + delete pe; + } +} +/*}}}*/ +/*FUNCTION Tria::DeleteResults {{{1*/ +void Tria::DeleteResults(void){ + + /*Delete and reinitialize results*/ + delete this->results; + this->results=new Results(); + +} +/*}}}*/ +/*FUNCTION Tria::GetNodeIndex {{{1*/ +int Tria::GetNodeIndex(Node* node){ + + ISSMASSERT(nodes); + for(int i=0;iGetParameterValue(&onbed,ElementOnBedEnum); + return onbed; +} +/*}}}*/ +/*FUNCTION Tria::IsOnShelf {{{1*/ +bool Tria::IsOnShelf(){ + + bool shelf; + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + return shelf; +} +/*}}}*/ +/*FUNCTION Tria::IsOnWater {{{1*/ +bool Tria::IsOnWater(){ + + bool water; + inputs->GetParameterValue(&water,ElementOnWaterEnum); + return water; +} +/*}}}*/ +/*FUNCTION Tria::GetSolutionFromInputs{{{1*/ +void Tria::GetSolutionFromInputs(Vec solution){ + + /*retrive parameters: */ + int analysis_type; + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */ + if (analysis_type==DiagnosticHorizAnalysisEnum) + GetSolutionFromInputsDiagnosticHoriz(solution); + else if (analysis_type==DiagnosticHutterAnalysisEnum) + GetSolutionFromInputsDiagnosticHutter(solution); + else + ISSMERROR("analysis: %s not supported yet",EnumToString(analysis_type)); + +} +/*}}}*/ +/*FUNCTION Tria::GetVectorFromInputs{{{1*/ +void Tria::GetVectorFromInputs(Vec vector,int NameEnum){ + + int doflist1[NUMVERTICES]; + + /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */ + for(int i=0;iinputs->Size();i++){ + Input* input=(Input*)this->inputs->GetObjectByOffset(i); + if(input->EnumType()==NameEnum){ + /*We found the enum. Use its values to fill into the vector, using the vertices ids: */ + this->GetDofList1(&doflist1[0]); + input->GetVectorFromInputs(vector,&doflist1[0]); + break; + } + } +} +/*}}}*/ +/*FUNCTION Tria::Gradj {{{1*/ +void Tria::Gradj(Vec gradient,int control_type){ + + /*If on water, grad = 0: */ + if(IsOnWater())return; + + switch(control_type){ + case DragCoefficientEnum: + GradjDrag(gradient); + break; + case RheologyBbarEnum: + GradjB(gradient); + break; + case DhDtEnum: + GradjDhDt(gradient); + break; + case VxEnum: + GradjVx(gradient); + break; + case VyEnum: + GradjVy(gradient); + break; + default: + ISSMERROR("%s%i","control type not supported yet: ",control_type); + } +} +/*}}}*/ +/*FUNCTION Tria::GradjB{{{1*/ +void Tria::GradjB(Vec gradient){ + + /*Intermediaries*/ + int i,ig; + int doflist[NUMVERTICES]; + double vx,vy,lambda,mu,thickness,Jdet; + double cm_noisedmp,viscosity_complement; + double dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2]; + double xyz_list[NUMVERTICES][3]; + double basis[3],epsilon[3]; + double dbasis[NDOF2][NUMVERTICES]; + double grad_g[NUMVERTICES]; + double grad[NUMVERTICES]={0.0}; + GaussTria *gauss = NULL; + + /*retrieve some parameters: */ + this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + GetDofList1(&doflist[0]); + + /*Retrieve all inputs*/ + Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); + Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); + Input* adjointx_input=inputs->GetInput(AdjointxEnum); ISSMASSERT(adjointx_input); + Input* adjointy_input=inputs->GetInput(AdjointyEnum); ISSMASSERT(adjointy_input); + Input* rheologyb_input=matice->inputs->GetInput(RheologyBbarEnum); ISSMASSERT(rheologyb_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(4); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + thickness_input->GetParameterValue(&thickness,gauss); + rheologyb_input->GetParameterDerivativeValue(&dB[0],&xyz_list[0][0],gauss); + vx_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss); + vy_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss); + adjointx_input->GetParameterDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss); + adjointy_input->GetParameterDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss); + + this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(basis,gauss); + GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss); + + /*standard gradient dJ/dki*/ + for (i=0;iweight*basis[i]; + } + /*Add regularization term*/ + for (i=0;iweight*(dbasis[0][i]*dB[0]+dbasis[1][i]*dB[1]); + for(i=0;iFindParam(&analysis_type,AnalysisTypeEnum); + + /*retrieve some parameters ands return if iceshelf: */ + this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); + if(IsOnShelf())return; + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + GetDofList1(&doflist1[0]); + + /*Build frictoin element, needed later: */ + inputs->GetParameterValue(&drag_type,DragTypeEnum); + friction=new Friction("2d",inputs,matpar,analysis_type); + + /*Retrieve all inputs we will be needing: */ + Input* adjointx_input=inputs->GetInput(AdjointxEnum); ISSMASSERT(adjointx_input); + Input* adjointy_input=inputs->GetInput(AdjointyEnum); ISSMASSERT(adjointy_input); + Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); + Input* dragcoefficient_input=inputs->GetInput(DragCoefficientEnum); ISSMASSERT(dragcoefficient_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(4); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(l1l2l3, gauss); + GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss); + + /*Build alpha_complement_list: */ + if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum); + else alpha_complement=0; + + dragcoefficient_input->GetParameterValue(&drag, gauss); + adjointx_input->GetParameterValue(&lambda, gauss); + adjointy_input->GetParameterValue(&mu, gauss); + vx_input->GetParameterValue(&vx,gauss); + vy_input->GetParameterValue(&vy,gauss); + dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for (i=0;iweight*l1l2l3[i]; + + //noise dampening d/dki(1/2*(dk/dx)^2) + grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); + } + + /*Add gradje_g_gaussian vector to gradje_g: */ + for( i=0; iGetInput(AdjointEnum); ISSMASSERT(adjoint_input); + Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(l1l2l3, gauss); + + adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); + thickness_input->GetParameterValue(&thickness, gauss); + + for(i=0;iweight*l1l2l3[i]; + } + + VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); + + /*Clean up and return*/ + delete gauss; +} +/*}}}*/ +/*FUNCTION Tria::GradjVy{{{1*/ +void Tria::GradjVy(Vec gradient){ + + /*Intermediaries*/ + int i,ig; + int doflist1[NUMVERTICES]; + double thickness,Jdet; + double l1l2l3[3]; + double Dlambda[2]; + double xyz_list[NUMVERTICES][3]; + double grade_g[NUMVERTICES] = {0.0}; + GaussTria *gauss = NULL; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + GetDofList1(&doflist1[0]); + + /*Retrieve all inputs we will be needing: */ + Input* adjoint_input=inputs->GetInput(AdjointEnum); ISSMASSERT(adjoint_input); + Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss); + thickness_input->GetParameterValue(&thickness, gauss); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctions(l1l2l3, gauss); + + for(i=0;iweight*l1l2l3[i]; + } + + VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); + + /*Clean up and return*/ + delete gauss; +} +/*}}}*/ +/*FUNCTION Tria::InputControlUpdate{{{1*/ +void Tria::InputControlUpdate(double scalar,bool save_parameter){ + + /*Intermediary*/ + int num_controls; + int* control_type=NULL; + Input* input=NULL; + double cm_min,cm_max; + + /*retrieve some parameters: */ + this->parameters->FindParam(&cm_min,CmMinEnum); + this->parameters->FindParam(&cm_max,CmMaxEnum); + this->parameters->FindParam(&num_controls,NumControlsEnum); + this->parameters->FindParam(&control_type,NULL,ControlTypeEnum); + + for(int i=0;iinputs->GetInput(control_type[i]); ISSMASSERT(input); + } + else{ + input=(Input*)this->inputs->GetInput(control_type[i]); ISSMASSERT(input); + } + + if (input->Enum()!=ControlInputEnum){ + ISSMERROR("input %s is not a ControlInput",EnumToString(control_type[i])); + } + + ((ControlInput*)input)->UpdateValue(scalar); + input->Constrain(cm_min,cm_max); + if (save_parameter) ((ControlInput*)input)->SaveValue(); + + } + + /*Clean up and return*/ + xfree((void**)&control_type); +} +/*}}}*/ +/*FUNCTION Tria::InputConvergence{{{1*/ +bool Tria::InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){ + + bool converged=true; + int i; + Input** new_inputs=NULL; + Input** old_inputs=NULL; + + new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs + old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs + + for(i=0;iinputs->GetInput(enums[2*i+0]); + old_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+1]); + if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumToString(enums[2*i+0])); + if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumToString(enums[2*i+0])); + } + + /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/ + for(i=0;icriterionvalues[i]) converged=false; + } + + /*clean up and return*/ + xfree((void**)&new_inputs); + xfree((void**)&old_inputs); + return converged; +} +/*}}}*/ +/*FUNCTION Tria::InputDepthAverageAtBase {{{1*/ +void Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){ + + /*New input*/ + Input* oldinput=NULL; + Input* newinput=NULL; + + /*copy input of enum_type*/ + if (object_enum==ElementsEnum) + oldinput=(Input*)this->inputs->GetInput(enum_type); + else if (object_enum==MaterialsEnum) + oldinput=(Input*)this->matice->inputs->GetInput(enum_type); + else + ISSMERROR("object %s not supported yet",EnumToString(object_enum)); + if(!oldinput)ISSMERROR("%s%s"," could not find old input with enum: ",EnumToString(enum_type)); + newinput=(Input*)oldinput->copy(); + + /*Assign new name (average)*/ + newinput->ChangeEnum(average_enum_type); + + /*Add new input to current element*/ + if (object_enum==ElementsEnum) + this->inputs->AddInput((Input*)newinput); + else if (object_enum==MaterialsEnum) + this->matice->inputs->AddInput((Input*)newinput); + else + ISSMERROR("object %s not supported yet",EnumToString(object_enum)); +} +/*}}}*/ +/*FUNCTION Tria::InputDuplicate{{{1*/ +void Tria::InputDuplicate(int original_enum,int new_enum){ + + /*Call inputs method*/ + if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum); + +} +/*}}}*/ +/*FUNCTION Tria::InputScale{{{1*/ +void Tria::InputScale(int enum_type,double scale_factor){ + + Input* input=NULL; + + /*Make a copy of the original input: */ + input=(Input*)this->inputs->GetInput(enum_type); + if(!input)ISSMERROR(" could not find old input with enum: %s",EnumToString(enum_type)); + + /*Scale: */ + input->Scale(scale_factor); +} +/*}}}*/ +/*FUNCTION Tria::InputArtificialNoise{{{1*/ +void Tria::InputArtificialNoise(int enum_type,double min,double max){ + + Input* input=NULL; + + /*Make a copy of the original input: */ + input=(Input*)this->inputs->GetInput(enum_type); + if(!input)ISSMERROR(" could not find old input with enum: %s",EnumToString(enum_type)); + + /*ArtificialNoise: */ + input->ArtificialNoise(min,max); +} +/*}}}*/ +/*FUNCTION Tria::InputToResult{{{1*/ +void Tria::InputToResult(int enum_type,int step,double time){ + + int i; + Input *input = NULL; + + /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ + if (enum_type==RheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type); + else input=this->inputs->GetInput(enum_type); + if (!input) ISSMERROR("Input %s not found in tria->inputs",EnumToString(enum_type)); + + /*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result + * object out of the input, with the additional step and time information: */ + this->results->AddObject((Object*)input->SpawnResult(step,time)); + if(input->Enum()==ControlInputEnum) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time)); +} +/*}}}*/ +/*FUNCTION Tria::MassFlux {{{1*/ +double Tria::MassFlux( double* segment,bool process_units){ + + const int numdofs=2; + + int i; + double mass_flux=0; + double xyz_list[NUMVERTICES][3]; + double normal[2]; + double length,rho_ice; + double x1,y1,x2,y2,h1,h2; + double vx1,vx2,vy1,vy2; + GaussTria* gauss_1=NULL; + GaussTria* gauss_2=NULL; + + /*Get material parameters :*/ + rho_ice=matpar->GetRhoIce(); + + /*First off, check that this segment belongs to this element: */ + if ((int)*(segment+4)!=this->id)ISSMERROR("%s%i%s%i","error message: segment with id ",(int)*(segment+4)," does not belong to element with id:",this->id); + + /*Recover segment node locations: */ + x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3); + + /*Get xyz list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + + /*get area coordinates of 0 and 1 locations: */ + gauss_1=new GaussTria(); + gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]); + gauss_2=new GaussTria(); + gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]); + + normal[0]=cos(atan2(x1-x2,y2-y1)); + normal[1]=sin(atan2(x1-x2,y2-y1)); + + length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2)); + + Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); + Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); + + thickness_input->GetParameterValue(&h1, gauss_1); + thickness_input->GetParameterValue(&h2, gauss_2); + vx_input->GetParameterValue(&vx1,gauss_1); + vx_input->GetParameterValue(&vx2,gauss_2); + vy_input->GetParameterValue(&vy1,gauss_1); + vy_input->GetParameterValue(&vy2,gauss_2); + + mass_flux= rho_ice*length*( + (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+ + (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1] + ); + + /*Process units: */ + mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum,this->parameters); + + /*clean up and return:*/ + delete gauss_1; + delete gauss_2; + return mass_flux; +} +/*}}}*/ +/*FUNCTION Tria::MaxAbsVx{{{1*/ +void Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){ + + /*Get maximum:*/ + double maxabsvx=this->inputs->MaxAbs(VxEnum); + + /*process units if requested: */ + if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum,this->parameters); + + /*Assign output pointers:*/ + *pmaxabsvx=maxabsvx; +} +/*}}}*/ +/*FUNCTION Tria::MaxAbsVy{{{1*/ +void Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){ + + /*Get maximum:*/ + double maxabsvy=this->inputs->MaxAbs(VyEnum); + + /*process units if requested: */ + if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum,this->parameters); + + /*Assign output pointers:*/ + *pmaxabsvy=maxabsvy; +} +/*}}}*/ +/*FUNCTION Tria::MaxAbsVz{{{1*/ +void Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){ + + /*Get maximum:*/ + double maxabsvz=this->inputs->MaxAbs(VzEnum); + + /*process units if requested: */ + if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum,this->parameters); + + /*Assign output pointers:*/ + *pmaxabsvz=maxabsvz; +} +/*}}}*/ +/*FUNCTION Tria::MaxVel{{{1*/ +void Tria::MaxVel(double* pmaxvel, bool process_units){ + + /*Get maximum:*/ + double maxvel=this->inputs->Max(VelEnum); + + /*process units if requested: */ + if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum,this->parameters); + + /*Assign output pointers:*/ + *pmaxvel=maxvel; +} +/*}}}*/ +/*FUNCTION Tria::MaxVx{{{1*/ +void Tria::MaxVx(double* pmaxvx, bool process_units){ + + /*Get maximum:*/ + double maxvx=this->inputs->Max(VxEnum); + + /*process units if requested: */ + if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum,this->parameters); + + /*Assign output pointers:*/ + *pmaxvx=maxvx; +} +/*}}}*/ +/*FUNCTION Tria::MaxVy{{{1*/ +void Tria::MaxVy(double* pmaxvy, bool process_units){ + + /*Get maximum:*/ + double maxvy=this->inputs->Max(VyEnum); + + /*process units if requested: */ + if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum,this->parameters); + + /*Assign output pointers:*/ + *pmaxvy=maxvy; + +} +/*}}}*/ +/*FUNCTION Tria::MaxVz{{{1*/ +void Tria::MaxVz(double* pmaxvz, bool process_units){ + + /*Get maximum:*/ + double maxvz=this->inputs->Max(VzEnum); + + /*process units if requested: */ + if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum,this->parameters); + + /*Assign output pointers:*/ + *pmaxvz=maxvz; +} +/*}}}*/ +/*FUNCTION Tria::MinVel{{{1*/ +void Tria::MinVel(double* pminvel, bool process_units){ + + /*Get minimum:*/ + double minvel=this->inputs->Min(VelEnum); + + /*process units if requested: */ + if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum,this->parameters); + + /*Assign output pointers:*/ + *pminvel=minvel; +} +/*}}}*/ +/*FUNCTION Tria::MinVx{{{1*/ +void Tria::MinVx(double* pminvx, bool process_units){ + + /*Get minimum:*/ + double minvx=this->inputs->Min(VxEnum); + + /*process units if requested: */ + if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum,this->parameters); + + /*Assign output pointers:*/ + *pminvx=minvx; +} +/*}}}*/ +/*FUNCTION Tria::MinVy{{{1*/ +void Tria::MinVy(double* pminvy, bool process_units){ + + /*Get minimum:*/ + double minvy=this->inputs->Min(VyEnum); + + /*process units if requested: */ + if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum,this->parameters); + + /*Assign output pointers:*/ + *pminvy=minvy; +} +/*}}}*/ +/*FUNCTION Tria::MinVz{{{1*/ +void Tria::MinVz(double* pminvz, bool process_units){ + + /*Get minimum:*/ + double minvz=this->inputs->Min(VzEnum); + + /*process units if requested: */ + if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum,this->parameters); + + /*Assign output pointers:*/ + *pminvz=minvz; +} +/*}}}*/ +/*FUNCTION Tria::TimeAdapt{{{1*/ +double Tria::TimeAdapt(void){ + + /*intermediary: */ + int i; + double C,dt; + double dx,dy; + double maxx,minx; + double maxy,miny; + double maxabsvx,maxabsvy; + double xyz_list[NUMVERTICES][3]; + + /*get CFL coefficient:*/ + this->parameters->FindParam(&C,CflCoefficientEnum); + + /*Get for Vx and Vy, the max of abs value: */ + this->MaxAbsVx(&maxabsvx,false); + this->MaxAbsVy(&maxabsvy,false); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], this->nodes, NUMVERTICES); + + minx=xyz_list[0][0]; + maxx=xyz_list[0][0]; + miny=xyz_list[0][1]; + maxy=xyz_list[0][1]; + + for(i=1;imaxx)maxx=xyz_list[i][0]; + if (xyz_list[i][1]maxy)maxy=xyz_list[i][1]; + } + dx=maxx-minx; + dy=maxy-miny; + + /*CFL criterion: */ + dt=C/(maxabsvy/dx+maxabsvy/dy); + + return dt; +} +/*}}}*/ +/*FUNCTION Tria::ThicknessAbsMisfit {{{1*/ +double Tria::ThicknessAbsMisfit(bool process_units){ + + /* Constants */ + const int numdof=1*NUMVERTICES; + + /*Intermediaries*/ + int i,ig; + double thickness,thicknessobs,weight; + double Jdet; + double Jelem = 0; + double xyz_list[NUMVERTICES][3]; + GaussTria *gauss = NULL; + + /*If on water, return 0: */ + if(IsOnWater())return 0; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + + /*Retrieve all inputs we will be needing: */ + Input* thickness_input =inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input); + Input* thicknessobs_input=inputs->GetInput(ThicknessObsEnum);ISSMASSERT(thicknessobs_input); + Input* weights_input =inputs->GetInput(WeightsEnum); ISSMASSERT(weights_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + + /*Get parameters at gauss point*/ + thickness_input->GetParameterValue(&thickness,gauss); + thicknessobs_input->GetParameterValue(&thicknessobs,gauss); + weights_input->GetParameterValue(&weight,gauss); + + /*compute ThicknessAbsMisfit*/ + Jelem+=0.5*pow(thickness-thicknessobs,2.0)*weight*Jdet*gauss->weight; + } + + /* clean up and Return: */ + delete gauss; + return Jelem; +} +/*}}}*/ +/*FUNCTION Tria::SurfaceAbsVelMisfit {{{1*/ +double Tria::SurfaceAbsVelMisfit(bool process_units){ + + const int numdof=NDOF2*NUMVERTICES; + + int i,ig; + double Jelem=0; + double velocity_mag,obs_velocity_mag; + double meanvel, epsvel,misfit,Jdet; + double xyz_list[NUMVERTICES][3]; + double vx_list[NUMVERTICES]; + double vy_list[NUMVERTICES]; + double obs_vx_list[NUMVERTICES]; + double obs_vy_list[NUMVERTICES]; + double misfit_square_list[NUMVERTICES]; + double misfit_list[NUMVERTICES]; + double weights_list[NUMVERTICES]; + GaussTria *gauss=NULL; + + /*If on water, return 0: */ + if(IsOnWater())return 0; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + + /* Recover input data: */ + GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); + GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); + GetParameterListOnVertices(&vx_list[0],VxEnum); + GetParameterListOnVertices(&vy_list[0],VyEnum); + GetParameterListOnVertices(&weights_list[0],WeightsEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&meanvel,MeanVelEnum); + this->parameters->FindParam(&epsvel,EpsVelEnum); + + /* Compute SurfaceAbsVelMisfit at the 3 nodes + * Here we integrate linearized functions: + * + * J(E) = int_E sum_{i=1}^3 J_i Phi_i + * + * where J_i are the misfits at the 3 nodes of the triangle + * Phi_i is the nodal function (P1) with respect to + * the vertex i + */ + + /*We are using an absolute misfit: + * + * 1 [ 2 2 ] + * J = --- | (u - u ) + (v - v ) | + * 2 [ obs obs ] + * + */ + for (i=0;iparameters); + + /*Apply weights to misfits*/ + for (i=0;ibegin();igend(); ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); + Jelem+=misfit*Jdet*gauss->weight; + } + + /*clean up and Return: */ + delete gauss; + return Jelem; +} +/*}}}*/ +/*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/ +double Tria::SurfaceRelVelMisfit(bool process_units){ + + const int numdof=2*NUMVERTICES; + + int i,ig; + double Jelem=0; + double scalex=1,scaley=1; + double meanvel, epsvel,misfit,Jdet; + double velocity_mag,obs_velocity_mag; + double xyz_list[NUMVERTICES][3]; + double vx_list[NUMVERTICES]; + double vy_list[NUMVERTICES]; + double obs_vx_list[NUMVERTICES]; + double obs_vy_list[NUMVERTICES]; + double misfit_square_list[NUMVERTICES]; + double misfit_list[NUMVERTICES]; + double weights_list[NUMVERTICES]; + GaussTria *gauss=NULL; + + /*If on water, return 0: */ + if(IsOnWater())return 0; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + + /* Recover input data: */ + GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); + GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); + GetParameterListOnVertices(&vx_list[0],VxEnum); + GetParameterListOnVertices(&vy_list[0],VyEnum); + GetParameterListOnVertices(&weights_list[0],WeightsEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&meanvel,MeanVelEnum); + this->parameters->FindParam(&epsvel,EpsVelEnum); + + /* Compute SurfaceRelVelMisfit at the 3 nodes + * Here we integrate linearized functions: + * + * J(E) = int_E sum_{i=1}^3 J_i Phi_i + * + * where J_i are the misfits at the 3 nodes of the triangle + * Phi_i is the nodal function (P1) with respect to + * the vertex i + */ + + /*We are using a relative misfit: + * + * 1 [ \bar{v}^2 2 \bar{v}^2 2 ] + * J = --- | ------------- (u - u ) + ------------- (v - v ) | + * 2 [ (u + eps)^2 obs (v + eps)^2 obs ] + * obs obs + */ + for (i=0;iparameters); + + /*Apply weights to misfits*/ + for (i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); + Jelem+=misfit*Jdet*gauss->weight; + } + + /*clean up and Return: */ + delete gauss; + return Jelem; +} +/*}}}*/ +/*FUNCTION Tria::SurfaceLogVelMisfit {{{1*/ +double Tria::SurfaceLogVelMisfit(bool process_units){ + + const int numdof=NDOF2*NUMVERTICES; + + int i,ig; + double Jelem=0; + double scalex=1,scaley=1; + double meanvel, epsvel,misfit,Jdet; + double velocity_mag,obs_velocity_mag; + double xyz_list[NUMVERTICES][3]; + double vx_list[NUMVERTICES]; + double vy_list[NUMVERTICES]; + double obs_vx_list[NUMVERTICES]; + double obs_vy_list[NUMVERTICES]; + double misfit_square_list[NUMVERTICES]; + double misfit_list[NUMVERTICES]; + double weights_list[NUMVERTICES]; + GaussTria *gauss=NULL; + + /*If on water, return 0: */ + if(IsOnWater())return 0; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + + /* Recover input data: */ + GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); + GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); + GetParameterListOnVertices(&vx_list[0],VxEnum); + GetParameterListOnVertices(&vy_list[0],VyEnum); + GetParameterListOnVertices(&weights_list[0],WeightsEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&meanvel,MeanVelEnum); + this->parameters->FindParam(&epsvel,EpsVelEnum); + + /* Compute SurfaceLogVelMisfit at the 3 nodes + * Here we integrate linearized functions: + * + * J(E) = int_E sum_{i=1}^3 J_i Phi_i + * + * where J_i are the misfits at the 3 nodes of the triangle + * Phi_i is the nodal function (P1) with respect to + * the vertex i + */ + + /*We are using a logarithmic misfit: + * + * [ vel + eps ] 2 + * J = 4 \bar{v}^2 | log ( ----------- ) | + * [ vel + eps ] + * obs + */ + for (i=0;iparameters); + + /*Apply weights to misfits*/ + for (i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); + Jelem+=misfit*Jdet*gauss->weight; + } + + /*clean-up and Return: */ + delete gauss; + return Jelem; +} +/*}}}*/ +/*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/ +double Tria::SurfaceLogVxVyMisfit(bool process_units){ + + const int numdof=NDOF2*NUMVERTICES; + + int i,ig; + int fit=-1; + double Jelem=0, S=0; + double scalex=1,scaley=1; + double meanvel, epsvel, misfit, Jdet; + double velocity_mag,obs_velocity_mag; + double xyz_list[NUMVERTICES][3]; + double vx_list[NUMVERTICES]; + double vy_list[NUMVERTICES]; + double obs_vx_list[NUMVERTICES]; + double obs_vy_list[NUMVERTICES]; + double misfit_square_list[NUMVERTICES]; + double misfit_list[NUMVERTICES]; + double weights_list[NUMVERTICES]; + GaussTria *gauss=NULL; + + /*If on water, return 0: */ + if(IsOnWater())return 0; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + + /* Recover input data: */ + GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); + GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); + GetParameterListOnVertices(&vx_list[0],VxEnum); + GetParameterListOnVertices(&vy_list[0],VyEnum); + GetParameterListOnVertices(&weights_list[0],WeightsEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&meanvel,MeanVelEnum); + this->parameters->FindParam(&epsvel,EpsVelEnum); + + /* Compute SurfaceLogVxVyMisfit at the 3 nodes + * Here we integrate linearized functions: + * + * J(E) = int_E sum_{i=1}^3 J_i Phi_i + * + * where J_i are the misfits at the 3 nodes of the triangle + * Phi_i is the nodal function (P1) with respect to + * the vertex i + */ + + /*We are using an logarithmic 2 misfit: + * + * 1 [ |u| + eps 2 |v| + eps 2 ] + * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) | + * 2 [ |u |+ eps |v |+ eps ] + * obs obs + */ + for (i=0;iparameters); + + /*Apply weights to misfits*/ + for (i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); + Jelem+=misfit*Jdet*gauss->weight; + } + + /*clean-up and Return: */ + delete gauss; + return Jelem; +} +/*}}}*/ +/*FUNCTION Tria::SurfaceAverageVelMisfit {{{1*/ +double Tria::SurfaceAverageVelMisfit(bool process_units){ + + const int numdof=2*NUMVERTICES; + + int i,ig; + int fit=-1; + double Jelem=0,S=0; + double scalex=1, scaley=1; + double meanvel, epsvel,Jdet; + double velocity_mag,obs_velocity_mag,misfit; + double xyz_list[NUMVERTICES][3]; + double vx_list[NUMVERTICES]; + double vy_list[NUMVERTICES]; + double obs_vx_list[NUMVERTICES]; + double obs_vy_list[NUMVERTICES]; + double misfit_square_list[NUMVERTICES]; + double misfit_list[NUMVERTICES]; + double weights_list[NUMVERTICES]; + GaussTria *gauss=NULL; + + /*If on water, return 0: */ + if(IsOnWater())return 0; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + + /* Recover input data: */ + inputs->GetParameterValue(&S,SurfaceAreaEnum); + GetParameterListOnVertices(&obs_vx_list[0],VxObsEnum); + GetParameterListOnVertices(&obs_vy_list[0],VyObsEnum); + GetParameterListOnVertices(&vx_list[0],VxEnum); + GetParameterListOnVertices(&vy_list[0],VyEnum); + GetParameterListOnVertices(&weights_list[0],WeightsEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&meanvel,MeanVelEnum); + this->parameters->FindParam(&epsvel,EpsVelEnum); + + /* Compute SurfaceAverageVelMisfit at the 3 nodes + * Here we integrate linearized functions: + * + * J(E) = int_E sum_{i=1}^3 J_i Phi_i + * + * where J_i are the misfits at the 3 nodes of the triangle + * Phi_i is the nodal function (P1) with respect to + * the vertex i + */ + + /*We are using a spacially average absolute misfit: + * + * 1 2 2 + * J = --- sqrt( (u - u ) + (v - v ) ) + * S obs obs + */ + for (i=0;iparameters); + + /*Take the square root, and scale by surface: */ + for (i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss); + Jelem+=misfit*Jdet*gauss->weight; + } + + /*clean-up and Return: */ + delete gauss; + return Jelem; +} +/*}}}*/ +/*FUNCTION Tria::PatchFill{{{1*/ +void Tria::PatchFill(int* prow, Patch* patch){ + + int i,row; + int vertices_ids[3]; + + /*recover pointer: */ + row=*prow; + + for(i=0;i<3;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch. + + for(i=0;iresults->Size();i++){ + ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); + + /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand + *it to the result object, to fill the rest: */ + patch->fillelementinfo(row,this->id,vertices_ids,3); + elementresult->PatchFill(row,patch); + + /*increment rower: */ + row++; + } + + /*Assign output pointers:*/ + *prow=row; +} +/*}}}*/ +/*FUNCTION Tria::PatchSize{{{1*/ +void Tria::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){ + + int i; + int numrows = 0; + int numnodes = 0; + + /*Go through all the results objects, and update the counters: */ + for (i=0;iresults->Size();i++){ + ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); + /*first, we have one more result: */ + numrows++; + /*now, how many vertices and how many nodal values for this result? :*/ + numnodes=elementresult->NumberOfNodalValues(); //ask result object. + } + + /*Assign output pointers:*/ + *pnumrows=numrows; + *pnumvertices=NUMVERTICES; + *pnumnodes=numnodes; +} +/*}}}*/ +/*FUNCTION Tria::ProcessResultsUnits{{{1*/ +void Tria::ProcessResultsUnits(void){ + + int i; + + for(i=0;iresults->Size();i++){ + ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); + elementresult->ProcessUnits(this->parameters); + } +} +/*}}}*/ +/*FUNCTION Tria::SurfaceArea {{{1*/ +double Tria::SurfaceArea(void){ + + int i; + double S; + double normal[3]; + double v13[3],v23[3]; + double xyz_list[NUMVERTICES][3]; + + /*If on water, return 0: */ + if(IsOnWater())return 0; + + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + + for (i=0;i<3;i++){ + v13[i]=xyz_list[0][i]-xyz_list[2][i]; + v23[i]=xyz_list[1][i]-xyz_list[2][i]; + } + + normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; + normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; + normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; + + S = 0.5 * sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2)); + + /*Return: */ + return S; +} +/*}}}*/ +/*FUNCTION Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){{{1*/ +void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index + + /*Intermediaries*/ + int i,j; + int tria_node_ids[3]; + int tria_vertex_ids[3]; + int tria_type; + double nodeinputs[3]; + + /*Checks if debuging*/ + /*{{{2*/ + ISSMASSERT(iomodel->elements); + /*}}}*/ + + /*Recover element type*/ + if ((analysis_type==PrognosticAnalysisEnum || analysis_type==BalancedthicknessAnalysisEnum) && iomodel->prognostic_DG){ + /*P1 Discontinuous Galerkin*/ + tria_type=P1DGEnum; + } + else{ + /*P1 Continuous Galerkin*/ + tria_type=P1Enum; + } + this->SetElementType(tria_type,analysis_counter); + + /*Recover vertices ids needed to initialize inputs*/ + for(i=0;i<3;i++){ + tria_vertex_ids[i]=(int)iomodel->elements[3*index+i]; //ids for vertices are in the elements array from Matlab + } + + /*Recover nodes ids needed to initialize the node hook.*/ + if (tria_type==P1DGEnum){ + /*Discontinuous Galerkin*/ + tria_node_ids[0]=iomodel->nodecounter+3*index+1; + tria_node_ids[1]=iomodel->nodecounter+3*index+2; + tria_node_ids[2]=iomodel->nodecounter+3*index+3; + } + else{ + /*Continuous Galerkin*/ + for(i=0;i<3;i++){ + tria_node_ids[i]=iomodel->nodecounter+(int)*(iomodel->elements+3*index+i); //ids for vertices are in the elements array from Matlab + } + } + + /*hooks: */ + this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type + + /*Fill with IoModel*/ + this->InputUpdateFromIoModel(index,iomodel); + + /*Defaults if not provided in iomodel*/ + switch(analysis_type){ + + case DiagnosticHorizAnalysisEnum: + + /*default vx,vy and vz: either observation or 0 */ + if(!iomodel->vx){ + if (iomodel->vx_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts; + else for(i=0;i<3;i++)nodeinputs[i]=0; + this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(VxOldEnum,nodeinputs)); + if(iomodel->qmu_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVxEnum,nodeinputs)); + } + if(!iomodel->vy){ + if (iomodel->vy_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts; + else for(i=0;i<3;i++)nodeinputs[i]=0; + this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(VyOldEnum,nodeinputs)); + if(iomodel->qmu_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVyEnum,nodeinputs)); + } + if(!iomodel->vz){ + if (iomodel->vz_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_vertex_ids[i]-1]/iomodel->yts; + else for(i=0;i<3;i++)nodeinputs[i]=0; + this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(VzOldEnum,nodeinputs)); + if(iomodel->qmu_analysis) this->inputs->AddInput(new TriaVertexInput(QmuVzEnum,nodeinputs)); + } + if(!iomodel->pressure){ + for(i=0;i<3;i++)nodeinputs[i]=0; + if(iomodel->qmu_analysis){ + this->inputs->AddInput(new TriaVertexInput(PressureEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(QmuPressureEnum,nodeinputs)); + } + } + break; + + default: + /*No update for other solution types*/ + break; + + } + + //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. + this->parameters=NULL; +} +/*}}}*/ /*FUNCTION Tria::UpdateGeometry{{{1*/ void Tria::UpdateGeometry(void){ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 6215) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 6216) @@ -66,4 +66,5 @@ void InputUpdateFromConstant(int constant, int name); void InputUpdateFromConstant(bool constant, int name); + void InputUpdateFromIoModel(int index, IoModel* iomodel); /*}}}*/ /*Element virtual functions definitions: {{{1*/ @@ -122,5 +123,4 @@ double SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); - void Update(int index, IoModel* iomodel); void UpdateGeometry(void); double TimeAdapt(); Index: /issm/trunk/src/c/objects/Loads/Icefront.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Icefront.h (revision 6215) +++ /issm/trunk/src/c/objects/Loads/Icefront.h (revision 6216) @@ -66,4 +66,5 @@ void InputUpdateFromConstant(bool constant, int name); void InputUpdateFromSolution(double* solution); + void InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");}; /*}}}*/ /*Load virtual functions definitions: {{{1*/ Index: /issm/trunk/src/c/objects/Loads/Numericalflux.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Numericalflux.h (revision 6215) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.h (revision 6216) @@ -62,4 +62,5 @@ void InputUpdateFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");} void InputUpdateFromSolution(double* solution){ISSMERROR("Not implemented yet!");} + void InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");}; /*}}}*/ /*Load virtual functions definitions: {{{1*/ Index: /issm/trunk/src/c/objects/Loads/Pengrid.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Pengrid.h (revision 6215) +++ /issm/trunk/src/c/objects/Loads/Pengrid.h (revision 6216) @@ -67,4 +67,5 @@ void InputUpdateFromConstant(bool constant, int name); void InputUpdateFromSolution(double* solution); + void InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");}; /*}}}*/ /*Load virtual functions definitions: {{{1*/ Index: /issm/trunk/src/c/objects/Loads/Penpair.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Penpair.h (revision 6215) +++ /issm/trunk/src/c/objects/Loads/Penpair.h (revision 6216) @@ -54,4 +54,5 @@ void InputUpdateFromConstant(bool constant, int name); void InputUpdateFromSolution(double* solution){ISSMERROR("Not implemented yet!");} + void InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");}; /*}}}*/ /*Load virtual functions definitions: {{{1*/ Index: /issm/trunk/src/c/objects/Loads/Riftfront.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Riftfront.h (revision 6215) +++ /issm/trunk/src/c/objects/Loads/Riftfront.h (revision 6216) @@ -73,4 +73,5 @@ void InputUpdateFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");} void InputUpdateFromSolution(double* solution){ISSMERROR("Not implemented yet!");} + void InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");}; /*}}}*/ /*Load virtual functions definitions: {{{1*/ Index: /issm/trunk/src/c/objects/Materials/Material.h =================================================================== --- /issm/trunk/src/c/objects/Materials/Material.h (revision 6215) +++ /issm/trunk/src/c/objects/Materials/Material.h (revision 6216) @@ -22,5 +22,4 @@ virtual void InputDuplicate(int original_enum,int new_enum)=0; virtual void Configure(Elements* elements)=0; - virtual void Update(int index, IoModel* iomodel)=0; }; Index: /issm/trunk/src/c/objects/Materials/Matice.cpp =================================================================== --- /issm/trunk/src/c/objects/Materials/Matice.cpp (revision 6215) +++ /issm/trunk/src/c/objects/Materials/Matice.cpp (revision 6216) @@ -248,59 +248,4 @@ return matice; -} -/*}}}*/ -/*FUNCTION Matice::Update{{{1*/ -void Matice::Update(int index, IoModel* iomodel){ - - int i,j; - - /*if 2d*/ - if(iomodel->dim==2){ - - /*Intermediaries*/ - const int num_vertices = 3; //Tria has 3 vertices - double nodeinputs[num_vertices]; - - /*Control Inputs*/ - if (iomodel->control_analysis && iomodel->control_type){ - for(i=0;inum_control_type;i++){ - switch((int)iomodel->control_type[i]){ - case RheologyBbarEnum: - if (iomodel->rheology_B){ - for(j=0;jrheology_B[int(iomodel->elements[num_vertices*index+j]-1)]; - this->inputs->AddInput(new ControlInput(RheologyBbarEnum,TriaVertexInputEnum,nodeinputs)); - } - break; - } - } - } - } - - /*if 3d*/ - else if(iomodel->dim==3){ - - /*Intermediaries*/ - const int num_vertices = 6; //Penta has 6 vertices - double nodeinputs[num_vertices]; - - /*Control Inputs*/ - if (iomodel->control_analysis && iomodel->control_type){ - for(i=0;inum_control_type;i++){ - switch((int)iomodel->control_type[i]){ - case RheologyBbarEnum: - if (iomodel->rheology_B){ - for(j=0;jrheology_B[int(iomodel->elements[num_vertices*index+j]-1)]; - this->inputs->AddInput(new ControlInput(RheologyBEnum,PentaVertexInputEnum,nodeinputs)); - } - break; - } - } - } - } - else{ - ISSMERROR(" Mesh type not supported yet!"); - } - - return; } /*}}}*/ @@ -721,4 +666,59 @@ } /*}}}*/ +/*FUNCTION Matice::InputUpdateFromIoModel{{{1*/ +void Matice::InputUpdateFromIoModel(int index, IoModel* iomodel){ + + int i,j; + + /*if 2d*/ + if(iomodel->dim==2){ + + /*Intermediaries*/ + const int num_vertices = 3; //Tria has 3 vertices + double nodeinputs[num_vertices]; + + /*Control Inputs*/ + if (iomodel->control_analysis && iomodel->control_type){ + for(i=0;inum_control_type;i++){ + switch((int)iomodel->control_type[i]){ + case RheologyBbarEnum: + if (iomodel->rheology_B){ + for(j=0;jrheology_B[int(iomodel->elements[num_vertices*index+j]-1)]; + this->inputs->AddInput(new ControlInput(RheologyBbarEnum,TriaVertexInputEnum,nodeinputs)); + } + break; + } + } + } + } + + /*if 3d*/ + else if(iomodel->dim==3){ + + /*Intermediaries*/ + const int num_vertices = 6; //Penta has 6 vertices + double nodeinputs[num_vertices]; + + /*Control Inputs*/ + if (iomodel->control_analysis && iomodel->control_type){ + for(i=0;inum_control_type;i++){ + switch((int)iomodel->control_type[i]){ + case RheologyBbarEnum: + if (iomodel->rheology_B){ + for(j=0;jrheology_B[int(iomodel->elements[num_vertices*index+j]-1)]; + this->inputs->AddInput(new ControlInput(RheologyBEnum,PentaVertexInputEnum,nodeinputs)); + } + break; + } + } + } + } + else{ + ISSMERROR(" Mesh type not supported yet!"); + } + + return; +} +/*}}}*/ /*FUNCTION Matice::IsInput{{{1*/ bool Matice::IsInput(int name){ Index: /issm/trunk/src/c/objects/Materials/Matice.h =================================================================== --- /issm/trunk/src/c/objects/Materials/Matice.h (revision 6215) +++ /issm/trunk/src/c/objects/Materials/Matice.h (revision 6216) @@ -52,5 +52,5 @@ void InputUpdateFromConstant(bool constant, int name); void InputUpdateFromSolution(double* solution); - void Update(int index, IoModel* iomodel); + void InputUpdateFromIoModel(int index, IoModel* iomodel); /*}}}*/ /*Material virtual functions resolution: {{{1*/ Index: /issm/trunk/src/c/objects/Materials/Matpar.h =================================================================== --- /issm/trunk/src/c/objects/Materials/Matpar.h (revision 6215) +++ /issm/trunk/src/c/objects/Materials/Matpar.h (revision 6216) @@ -56,5 +56,5 @@ void InputUpdateFromConstant(bool constant, int name); void InputUpdateFromSolution(double* solution); - void Update(int index, IoModel* iomodel){ISSMERROR("not implemented yet");}; + void InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("not implemented yet");}; /*}}}*/ /*Material virtual functions resolution: {{{1*/ Index: /issm/trunk/src/c/objects/Node.h =================================================================== --- /issm/trunk/src/c/objects/Node.h (revision 6215) +++ /issm/trunk/src/c/objects/Node.h (revision 6216) @@ -62,4 +62,5 @@ void InputUpdateFromConstant(bool constant, int name); void InputUpdateFromSolution(double* solution){ISSMERROR("Not implemented yet!");} + void InputUpdateFromIoModel(int index, IoModel* iomodel){ISSMERROR("Not implemented yet!");} /*}}}*/ /*Node numerical routines {{{1*/ Index: /issm/trunk/src/c/objects/Update.h =================================================================== --- /issm/trunk/src/c/objects/Update.h (revision 6215) +++ /issm/trunk/src/c/objects/Update.h (revision 6216) @@ -25,4 +25,5 @@ virtual void InputUpdateFromConstant(bool constant, int name)=0; virtual void InputUpdateFromSolution(double* solution)=0; + virtual void InputUpdateFromIoModel(int index, IoModel* iomodel)=0; };