Index: /issm/trunk/src/c/objects/Elements/Beam.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Beam.cpp (revision 4249) +++ /issm/trunk/src/c/objects/Elements/Beam.cpp (revision 4250) @@ -31,11 +31,34 @@ } /*}}}*/ -/*FUNCTION void Beam::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);{{{1*/ -void Beam::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ - -/*Object virtual functions definitions: */ + +/*Object virtual functions definitions:*/ +/*FUNCTION Beam::copy{{{1*/ +Object* Beam::copy() { + + int i; + Beam* beam=NULL; + + beam=new Beam(); + + /*copy fields: */ + beam->id=this->id; + if(this->inputs){ + beam->inputs=(Inputs*)this->inputs->Copy(); + } + else{ + beam->inputs=new Inputs(); + } + /*point parameters: */ + beam->parameters=this->parameters; + + /*pointers: */ + beam->nodes=(Node**)xmalloc(2*sizeof(Node*)); + for(i=0;i<2;i++)beam->nodes[i]=this->nodes[i]; + beam->matice=this->matice; + beam->matpar=this->matpar; + + return beam; +} +/*}}}*/ /*FUNCTION Beam::DeepEcho{{{1*/ void Beam::DeepEcho(void){ @@ -58,70 +81,7 @@ } /*}}}*/ -/*FUNCTION Beam::Id{{{1*/ -int Beam::Id(void){ return id; } -/*}}}*/ -/*FUNCTION Beam::MyRank{{{1*/ -int Beam::MyRank(void){ - extern int my_rank; - return my_rank; -} -/*}}}*/ -/*FUNCTION Beam::Marshall{{{1*/ -void Beam::Marshall(char** pmarshalled_dataset){ - ISSMERROR("not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::MarshallSize{{{1*/ -int Beam::MarshallSize(){ - ISSMERROR("not supported yet!"); -} -/*}}}*/ /*FUNCTION Beam::Demarshall{{{1*/ void Beam::Demarshall(char** pmarshalled_dataset){ ISSMERROR("not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::Enum{{{1*/ -int Beam::Enum(void){ - - return BeamEnum; - -} -/*}}}*/ -/*FUNCTION Beam::copy{{{1*/ -Object* Beam::copy() { - - int i; - Beam* beam=NULL; - - beam=new Beam(); - - /*copy fields: */ - beam->id=this->id; - if(this->inputs){ - beam->inputs=(Inputs*)this->inputs->Copy(); - } - else{ - beam->inputs=new Inputs(); - } - /*point parameters: */ - beam->parameters=this->parameters; - - /*pointers: */ - beam->nodes=(Node**)xmalloc(2*sizeof(Node*)); - for(i=0;i<2;i++)beam->nodes[i]=this->nodes[i]; - beam->matice=this->matice; - beam->matpar=this->matpar; - - return beam; -} -/*}}}*/ - -/*Beam management*/ -/*FUNCTION Beam::Configure {{{1*/ -void Beam::Configure(Elements* elementsin,Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ - - ISSMERROR(" not supported yet!"); - } /*}}}*/ @@ -146,11 +106,67 @@ } /*}}}*/ -/*FUNCTION Beam::IsInput{{{1*/ -bool Beam::IsInput(int name){ - if (name==SurfaceSlopeXEnum || - name==SurfaceSlopeYEnum){ - return true; - } - else return false; +/*FUNCTION Beam::Enum{{{1*/ +int Beam::Enum(void){ + + return BeamEnum; + +} +/*}}}*/ +/*FUNCTION Beam::Id{{{1*/ +int Beam::Id(void){ return id; } +/*}}}*/ +/*FUNCTION Beam::Marshall{{{1*/ +void Beam::Marshall(char** pmarshalled_dataset){ + ISSMERROR("not supported yet!"); +} +/*}}}*/ +/*FUNCTION Beam::MarshallSize{{{1*/ +int Beam::MarshallSize(){ + ISSMERROR("not supported yet!"); +} +/*}}}*/ +/*FUNCTION Beam::MyRank{{{1*/ +int Beam::MyRank(void){ + extern int my_rank; + return my_rank; +} +/*}}}*/ + +/*Update virtual functions definitions:*/ +/*FUNCTION Beam::InputUpdateFromVector(double* vector, int name, int type);{{{1*/ +void Beam::InputUpdateFromVector(double* vector, int name, int type){ + + /*Check that name is an element input*/ + if (!IsInput(name)) return; + switch(type){ + + case VertexEnum: + + /*New PentaVertexInpu*/ + double values[2]; + + /*Get values on the 6 vertices*/ + for (int i=0;i<2;i++){ + values[i]=vector[nodes[i]->GetVertexDof()]; + } + + /*update input*/ + this->inputs->AddInput(new BeamVertexInput(name,values)); + return; + + default: + + ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type)); + } +} +/*}}}*/ +/*FUNCTION Beam::InputUpdateFromVector(int* vector, int name, int type);{{{1*/ +void Beam::InputUpdateFromVector(int* vector, int name, int type){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Beam::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/ +void Beam::InputUpdateFromVector(bool* vector, int name, int type){ + ISSMERROR(" not supported yet!"); } /*}}}*/ @@ -160,21 +176,6 @@ } /*}}}*/ -/*FUNCTION Beam::GetSolutionFromInputs(Vec solution);{{{1*/ -void Beam::GetSolutionFromInputs(Vec solution){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::InputToResult(int enum_type,int step,double time){{{1*/ -void Beam::InputToResult(int enum_type,int step,double time){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::ProcessResultsUnits(void){{{1*/ -void Beam::ProcessResultsUnits(void){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ - -/*Object functions*/ + +/*Element virtual functions definitions:*/ /*FUNCTION Beam::ComputeBasalStress{{{1*/ void Beam::ComputeBasalStress(Vec eps){ @@ -232,4 +233,11 @@ } /*}}}*/ +/*FUNCTION Beam::Configure {{{1*/ +void Beam::Configure(Elements* elementsin,Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ + + ISSMERROR(" not supported yet!"); + +} +/*}}}*/ /*FUNCTION Beam::CostFunction{{{1*/ double Beam::CostFunction(void){ @@ -255,4 +263,563 @@ } /*}}}*/ +/*FUNCTION Beam::CreatePVector{{{1*/ +void Beam::CreatePVector(Vec pg){ + + int analysis_type,sub_analysis_type; + + /*retrive parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum); + + /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ + if (analysis_type==DiagnosticHutterAnalysisEnum) { + CreatePVectorDiagnosticHutter( pg); + } + else{ + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type)); + } + +} +/*}}}*/ +/*FUNCTION Beam::Du{{{1*/ +void Beam::Du(Vec){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Beam::GetBedList{{{1*/ +void Beam::GetBedList(double*){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Beam::GetMatPar{{{1*/ +void* Beam::GetMatPar(){ + + return matpar; +} +/*}}}*/ +/*FUNCTION Beam::GetNodes{{{1*/ +void Beam::GetNodes(void** vpnodes){ + int i; + Node** pnodes=NULL; + + /*recover nodes: */ + pnodes=(Node**)vpnodes; + + for(i=0;i<3;i++){ + pnodes[i]=nodes[i]; + } +} +/*}}}*/ +/*FUNCTION Beam::GetOnBed{{{1*/ +bool Beam::GetOnBed(){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Beam::GetShelf{{{1*/ +bool Beam::GetShelf(){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Beam::GetSolutionFromInputs(Vec solution);{{{1*/ +void Beam::GetSolutionFromInputs(Vec solution){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Beam::GetThicknessList{{{1*/ +void Beam::GetThicknessList(double* thickness_list){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Beam::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/ +void Beam::GetVectorFromInputs(Vec vector,int NameEnum){ + + int i; + const int numvertices=2; + int doflist1[numvertices]; + + /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */ + for(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 Beam::Gradj{{{1*/ +void Beam::Gradj(Vec gradient,int control_type){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Beam::GradjB{{{1*/ +void Beam::GradjB(Vec gradient){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Beam::GradjDrag{{{1*/ +void Beam::GradjDrag(Vec gradient){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Beam::InputAXPY(int YEnum, double scalar, int XEnum);{{{1*/ +void Beam::InputAXPY(int YEnum, double scalar, int XEnum){ + + Input* xinput=NULL; + Input* yinput=NULL; + + /*Find x and y inputs: */ + xinput=(Input*)this->inputs->GetInput(XEnum); + yinput=(Input*)this->inputs->GetInput(YEnum); + + /*some checks: */ + if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!"); + if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!"); + + /*Scale: */ + yinput->AXPY(xinput,scalar); +} +/*}}}*/ +/*FUNCTION Beam::InputControlConstrain(int control_type, double cm_min, double cm_max){{{1*/ +void Beam::InputControlConstrain(int control_type, double cm_min, double cm_max){ + + Input* input=NULL; + + /*Find input: */ + input=(Input*)this->inputs->GetInput(control_type); + + /*Do nothing if we don't find it: */ + if(!input)return; + + /*Constrain input using cm_min and cm_max: */ + input->Constrain(cm_min,cm_max); + +} +/*}}}*/ +/*FUNCTION Beam::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/ +void Beam::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){ + + int i; + Input** new_inputs=NULL; + Input** old_inputs=NULL; + int converged=1; + + 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 ",EnumAsString(enums[2*i+0])); + if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(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=0; + } + + /*Assign output pointers:*/ + *pconverged=converged; + +} +/*}}}*/ +/*FUNCTION Beam::InputDuplicate(int original_enum,int new_enum){{{1*/ +void Beam::InputDuplicate(int original_enum,int new_enum){ + + Input* original=NULL; + Input* copy=NULL; + + /*Make a copy of the original input: */ + original=(Input*)this->inputs->GetInput(original_enum); + copy=(Input*)original->copy(); + + /*Change copy enum to reinitialized_enum: */ + copy->ChangeEnum(new_enum); + + /*Add copy into inputs, it will wipe off the one already there: */ + inputs->AddObject((Input*)copy); +} +/*}}}*/ +/*FUNCTION Beam::InputScale(int enum_type,double scale_factor){{{1*/ +void Beam::InputScale(int enum_type,double scale_factor){ + + Input* input=NULL; + + /*Make a copy of the original input: */ + input=(Input*)this->inputs->GetInput(enum_type); + + /*Scale: */ + input->Scale(scale_factor); +} +/*}}}*/ +/*FUNCTION Beam::InputToResult(int enum_type,int step,double time){{{1*/ +void Beam::InputToResult(int enum_type,int step,double time){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Beam::MassFlux{{{1*/ +double Beam::MassFlux( double* segment){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Beam::MaxAbsVx(double* pmaxabsvx, bool process_units);{{{1*/ +void Beam::MaxAbsVx(double* pmaxabsvx, bool process_units){ + + int i; + int dim; + const int numgrids=2; + double gaussgrids[numgrids][2]={{0,1},{1,0}}; + double vx_values[numgrids]; + double maxabsvx; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); + + /*now, compute maximum:*/ + maxabsvx=fabs(vx_values[0]); + for(i=1;imaxabsvx)maxabsvx=fabs(vx_values[i]); + } + + /*Assign output pointers:*/ + *pmaxabsvx=maxabsvx; +} +/*}}}*/ +/*FUNCTION Beam::MaxAbsVy(double* pmaxabsvy, bool process_units);{{{1*/ +void Beam::MaxAbsVy(double* pmaxabsvy, bool process_units){ + + int i; + int dim; + const int numgrids=2; + double gaussgrids[numgrids][2]={{0,1},{1,0}}; + double vy_values[numgrids]; + double maxabsvy; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); + + /*now, compute maximum:*/ + maxabsvy=fabs(vy_values[0]); + for(i=1;imaxabsvy)maxabsvy=fabs(vy_values[i]); + } + + /*Assign output pointers:*/ + *pmaxabsvy=maxabsvy; +} +/*}}}*/ +/*FUNCTION Beam::MaxAbsVz(double* pmaxabsvz, bool process_units);{{{1*/ +void Beam::MaxAbsVz(double* pmaxabsvz, bool process_units){ + + int i; + int dim; + const int numgrids=2; + double gaussgrids[numgrids][2]={{0,1},{1,0}}; + double vz_values[numgrids]; + double maxabsvz; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); + + /*now, compute maximum:*/ + maxabsvz=fabs(vz_values[0]); + for(i=1;imaxabsvz)maxabsvz=fabs(vz_values[i]); + } + + /*Assign output pointers:*/ + *pmaxabsvz=maxabsvz; +} +/*}}}*/ +/*FUNCTION Beam::MaxVel(double* pmaxvel, bool process_units);{{{1*/ +void Beam::MaxVel(double* pmaxvel, bool process_units){ + + int i; + int dim; + const int numgrids=2; + double gaussgrids[numgrids][2]={{0,1},{1,0}}; + double vx_values[numgrids]; + double vy_values[numgrids]; + double vz_values[numgrids]; + double vel_values[numgrids]; + double maxvel; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); + inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); + if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); + + /*now, compute maximum of velocity :*/ + if(dim==2){ + for(i=0;imaxvel)maxvel=vel_values[i]; + } + + /*Assign output pointers:*/ + *pmaxvel=maxvel; + +} +/*}}}*/ +/*FUNCTION Beam::MaxVx(double* pmaxvx, bool process_units);{{{1*/ +void Beam::MaxVx(double* pmaxvx, bool process_units){ + + int i; + int dim; + const int numgrids=2; + double gaussgrids[numgrids][2]={{0,1},{1,0}}; + double vx_values[numgrids]; + double maxvx; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); + + /*now, compute maximum:*/ + maxvx=vx_values[0]; + for(i=1;imaxvx)maxvx=vx_values[i]; + } + + /*Assign output pointers:*/ + *pmaxvx=maxvx; + +} +/*}}}*/ +/*FUNCTION Beam::MaxVy(double* pmaxvy, bool process_units);{{{1*/ +void Beam::MaxVy(double* pmaxvy, bool process_units){ + + int i; + int dim; + const int numgrids=2; + double gaussgrids[numgrids][2]={{0,1},{1,0}}; + double vy_values[numgrids]; + double maxvy; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); + + /*now, compute maximum:*/ + maxvy=vy_values[0]; + for(i=1;imaxvy)maxvy=vy_values[i]; + } + + /*Assign output pointers:*/ + *pmaxvy=maxvy; + +} +/*}}}*/ +/*FUNCTION Beam::MaxVz(double* pmaxvz, bool process_units);{{{1*/ +void Beam::MaxVz(double* pmaxvz, bool process_units){ + + int i; + int dim; + const int numgrids=2; + double gaussgrids[numgrids][2]={{0,1},{1,0}}; + double vz_values[numgrids]; + double maxvz; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); + + /*now, compute maximum:*/ + maxvz=vz_values[0]; + for(i=1;imaxvz)maxvz=vz_values[i]; + } + + /*Assign output pointers:*/ + *pmaxvz=maxvz; + +} +/*}}}*/ +/*FUNCTION Beam::MinVel(double* pminvel, bool process_units);{{{1*/ +void Beam::MinVel(double* pminvel, bool process_units){ + + int i; + int dim; + const int numgrids=2; + double gaussgrids[numgrids][2]={{0,1},{1,0}}; + double vx_values[numgrids]; + double vy_values[numgrids]; + double vz_values[numgrids]; + double vel_values[numgrids]; + double minvel; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); + inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); + if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); + + /*now, compute minimum of velocity :*/ + if(dim==2){ + for(i=0;iFindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); + + /*now, compute minimum:*/ + minvx=vx_values[0]; + for(i=1;iFindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); + + /*now, compute minimum:*/ + minvy=vy_values[0]; + for(i=1;iFindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); + + /*now, compute minimum:*/ + minvz=vz_values[0]; + for(i=1;iFindParam(&analysis_type,AnalysisTypeEnum); - parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum); - - /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ - if (analysis_type==DiagnosticHutterAnalysisEnum) { - CreatePVectorDiagnosticHutter( pg); - } - else{ - ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type)); - } } @@ -455,14 +1003,4 @@ } /*}}}*/ -/*FUNCTION Beam::Du{{{1*/ -void Beam::Du(Vec){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::GetBedList{{{1*/ -void Beam::GetBedList(double*){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ /*FUNCTION Beam::GetDofList{{{1*/ void Beam::GetDofList(int* doflist,int* pnumberofdofspernode){ @@ -509,10 +1047,4 @@ } /*}}}*/ -/*FUNCTION Beam::GetMatPar{{{1*/ -void* Beam::GetMatPar(){ - - return matpar; -} -/*}}}*/ /*FUNCTION Beam::GetNodalFunctions{{{1*/ void Beam::GetNodalFunctions(double* l1l2, double gauss_coord){ @@ -527,22 +1059,4 @@ } /*}}}*/ -/*FUNCTION Beam::GetNodes{{{1*/ -void Beam::GetNodes(void** vpnodes){ - int i; - Node** pnodes=NULL; - - /*recover nodes: */ - pnodes=(Node**)vpnodes; - - for(i=0;i<3;i++){ - pnodes[i]=nodes[i]; - } -} -/*}}}*/ -/*FUNCTION Beam::GetOnBed{{{1*/ -bool Beam::GetOnBed(){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ /*FUNCTION Beam::GetParameterValue{{{1*/ void Beam::GetParameterValue(double* pvalue, double* value_list,double gauss_coord){ @@ -555,42 +1069,11 @@ } /*}}}*/ -/*FUNCTION Beam::GetShelf{{{1*/ -bool Beam::GetShelf(){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::GetThicknessList{{{1*/ -void Beam::GetThicknessList(double* thickness_list){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::Gradj{{{1*/ -void Beam::Gradj(Vec gradient,int control_type){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::GradjB{{{1*/ -void Beam::GradjB(Vec gradient){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::GradjDrag{{{1*/ -void Beam::GradjDrag(Vec gradient){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::MassFlux{{{1*/ -double Beam::MassFlux( double* segment){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::Misfit{{{1*/ -double Beam::Misfit(void){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::SurfaceArea{{{1*/ -double Beam::SurfaceArea(void){ - ISSMERROR(" not supported yet!"); +/*FUNCTION Beam::IsInput{{{1*/ +bool Beam::IsInput(int name){ + if (name==SurfaceSlopeXEnum || + name==SurfaceSlopeYEnum){ + return true; + } + else return false; } /*}}}*/ @@ -601,483 +1084,2 @@ } /*}}}1*/ -/*FUNCTION Beam::InputUpdateFromVector(double* vector, int name, int type);{{{1*/ -void Beam::InputUpdateFromVector(double* vector, int name, int type){ - - /*Check that name is an element input*/ - if (!IsInput(name)) return; - switch(type){ - - case VertexEnum: - - /*New PentaVertexInpu*/ - double values[2]; - - /*Get values on the 6 vertices*/ - for (int i=0;i<2;i++){ - values[i]=vector[nodes[i]->GetVertexDof()]; - } - - /*update input*/ - this->inputs->AddInput(new BeamVertexInput(name,values)); - return; - - default: - - ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type)); - } -} -/*}}}*/ -/*FUNCTION Beam::InputUpdateFromVector(int* vector, int name, int type);{{{1*/ -void Beam::InputUpdateFromVector(int* vector, int name, int type){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/ -void Beam::InputUpdateFromVector(bool* vector, int name, int type){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){{{1*/ -void Beam::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){ - - ISSMERROR(" not supported yet!"); - -} -/*}}}*/ -/*FUNCTION Beam::PatchFill(int* pcount, Patch* patch){{{1*/ -void Beam::PatchFill(int* pcount, Patch* patch){ - - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Beam::MinVel(double* pminvel, bool process_units);{{{1*/ -void Beam::MinVel(double* pminvel, bool process_units){ - - int i; - int dim; - const int numgrids=2; - double gaussgrids[numgrids][2]={{0,1},{1,0}}; - double vx_values[numgrids]; - double vy_values[numgrids]; - double vz_values[numgrids]; - double vel_values[numgrids]; - double minvel; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); - inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); - if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); - - /*now, compute minimum of velocity :*/ - if(dim==2){ - for(i=0;iFindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); - inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); - if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); - - /*now, compute maximum of velocity :*/ - if(dim==2){ - for(i=0;imaxvel)maxvel=vel_values[i]; - } - - /*Assign output pointers:*/ - *pmaxvel=maxvel; - -} -/*}}}*/ -/*FUNCTION Beam::MinVx(double* pminvx, bool process_units);{{{1*/ -void Beam::MinVx(double* pminvx, bool process_units){ - - int i; - int dim; - const int numgrids=2; - double gaussgrids[numgrids][2]={{0,1},{1,0}}; - double vx_values[numgrids]; - double minvx; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); - - /*now, compute minimum:*/ - minvx=vx_values[0]; - for(i=1;iFindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); - - /*now, compute maximum:*/ - maxvx=vx_values[0]; - for(i=1;imaxvx)maxvx=vx_values[i]; - } - - /*Assign output pointers:*/ - *pmaxvx=maxvx; - -} -/*}}}*/ -/*FUNCTION Beam::MaxAbsVx(double* pmaxabsvx, bool process_units);{{{1*/ -void Beam::MaxAbsVx(double* pmaxabsvx, bool process_units){ - - int i; - int dim; - const int numgrids=2; - double gaussgrids[numgrids][2]={{0,1},{1,0}}; - double vx_values[numgrids]; - double maxabsvx; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); - - /*now, compute maximum:*/ - maxabsvx=fabs(vx_values[0]); - for(i=1;imaxabsvx)maxabsvx=fabs(vx_values[i]); - } - - /*Assign output pointers:*/ - *pmaxabsvx=maxabsvx; -} -/*}}}*/ -/*FUNCTION Beam::MinVy(double* pminvy, bool process_units);{{{1*/ -void Beam::MinVy(double* pminvy, bool process_units){ - - int i; - int dim; - const int numgrids=2; - double gaussgrids[numgrids][2]={{0,1},{1,0}}; - double vy_values[numgrids]; - double minvy; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); - - /*now, compute minimum:*/ - minvy=vy_values[0]; - for(i=1;iFindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); - - /*now, compute maximum:*/ - maxvy=vy_values[0]; - for(i=1;imaxvy)maxvy=vy_values[i]; - } - - /*Assign output pointers:*/ - *pmaxvy=maxvy; - -} -/*}}}*/ -/*FUNCTION Beam::MaxAbsVy(double* pmaxabsvy, bool process_units);{{{1*/ -void Beam::MaxAbsVy(double* pmaxabsvy, bool process_units){ - - int i; - int dim; - const int numgrids=2; - double gaussgrids[numgrids][2]={{0,1},{1,0}}; - double vy_values[numgrids]; - double maxabsvy; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); - - /*now, compute maximum:*/ - maxabsvy=fabs(vy_values[0]); - for(i=1;imaxabsvy)maxabsvy=fabs(vy_values[i]); - } - - /*Assign output pointers:*/ - *pmaxabsvy=maxabsvy; -} -/*}}}*/ -/*FUNCTION Beam::MinVz(double* pminvz, bool process_units);{{{1*/ -void Beam::MinVz(double* pminvz, bool process_units){ - - int i; - int dim; - const int numgrids=2; - double gaussgrids[numgrids][2]={{0,1},{1,0}}; - double vz_values[numgrids]; - double minvz; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); - - /*now, compute minimum:*/ - minvz=vz_values[0]; - for(i=1;iFindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); - - /*now, compute maximum:*/ - maxvz=vz_values[0]; - for(i=1;imaxvz)maxvz=vz_values[i]; - } - - /*Assign output pointers:*/ - *pmaxvz=maxvz; - -} -/*}}}*/ -/*FUNCTION Beam::MaxAbsVz(double* pmaxabsvz, bool process_units);{{{1*/ -void Beam::MaxAbsVz(double* pmaxabsvz, bool process_units){ - - int i; - int dim; - const int numgrids=2; - double gaussgrids[numgrids][2]={{0,1},{1,0}}; - double vz_values[numgrids]; - double maxabsvz; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); - - /*now, compute maximum:*/ - maxabsvz=fabs(vz_values[0]); - for(i=1;imaxabsvz)maxabsvz=fabs(vz_values[i]); - } - - /*Assign output pointers:*/ - *pmaxabsvz=maxabsvz; -} -/*}}}*/ -/*FUNCTION Beam::InputDuplicate(int original_enum,int new_enum){{{1*/ -void Beam::InputDuplicate(int original_enum,int new_enum){ - - Input* original=NULL; - Input* copy=NULL; - - /*Make a copy of the original input: */ - original=(Input*)this->inputs->GetInput(original_enum); - copy=(Input*)original->copy(); - - /*Change copy enum to reinitialized_enum: */ - copy->ChangeEnum(new_enum); - - /*Add copy into inputs, it will wipe off the one already there: */ - inputs->AddObject((Input*)copy); -} -/*}}}*/ -/*FUNCTION Beam::InputScale(int enum_type,double scale_factor){{{1*/ -void Beam::InputScale(int enum_type,double scale_factor){ - - Input* input=NULL; - - /*Make a copy of the original input: */ - input=(Input*)this->inputs->GetInput(enum_type); - - /*Scale: */ - input->Scale(scale_factor); -} -/*}}}*/ -/*FUNCTION Beam::InputAXPY(int YEnum, double scalar, int XEnum);{{{1*/ -void Beam::InputAXPY(int YEnum, double scalar, int XEnum){ - - Input* xinput=NULL; - Input* yinput=NULL; - - /*Find x and y inputs: */ - xinput=(Input*)this->inputs->GetInput(XEnum); - yinput=(Input*)this->inputs->GetInput(YEnum); - - /*some checks: */ - if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!"); - if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!"); - - /*Scale: */ - yinput->AXPY(xinput,scalar); -} -/*}}}*/ -/*FUNCTION Beam::InputControlConstrain(int control_type, double cm_min, double cm_max){{{1*/ -void Beam::InputControlConstrain(int control_type, double cm_min, double cm_max){ - - Input* input=NULL; - - /*Find input: */ - input=(Input*)this->inputs->GetInput(control_type); - - /*Do nothing if we don't find it: */ - if(!input)return; - - /*Constrain input using cm_min and cm_max: */ - input->Constrain(cm_min,cm_max); - -} -/*}}}*/ -/*FUNCTION Beam::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/ -void Beam::GetVectorFromInputs(Vec vector,int NameEnum){ - - int i; - const int numvertices=2; - int doflist1[numvertices]; - - /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */ - for(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 Beam::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/ -void Beam::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){ - - int i; - Input** new_inputs=NULL; - Input** old_inputs=NULL; - int converged=1; - - 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 ",EnumAsString(enums[2*i+0])); - if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(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=0; - } - - /*Assign output pointers:*/ - *pconverged=converged; - -} -/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Beam.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Beam.h (revision 4249) +++ /issm/trunk/src/c/objects/Elements/Beam.h (revision 4250) @@ -42,13 +42,13 @@ /*}}}*/ /*Object virtual functions definitions: {{{1*/ - void Echo(); - void DeepEcho(); - int Id(); - int MyRank(); - void Marshall(char** pmarshalled_dataset); - int MarshallSize(); - void Demarshall(char** pmarshalled_dataset); - int Enum(); Object* copy(); + void DeepEcho(); + void Demarshall(char** pmarshalled_dataset); + void Echo(); + int Enum(); + int Id(); + void Marshall(char** pmarshalled_dataset); + int MarshallSize(); + int MyRank(); /*}}}*/ /*Update virtual functions resolution: {{{1*/ Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4249) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4250) @@ -88,151 +88,47 @@ } /*}}}*/ -/*FUNCTION Penta::Update(IoModel* iomodel,int analysis_counter,int analysis_type) {{{1*/ -void Penta::Update(int index,IoModel* iomodel,int analysis_counter,int analysis_type){ - - /*Intermediaries*/ - IssmInt i; - int penta_node_ids[6]; - 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 - } - - /*Recover nodes ids needed to initialize the node hook.*/ - for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook. - penta_node_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab - } - - /*hooks: */ - this->SetHookNodes(penta_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type - - //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]; - 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->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->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->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)); - } - - /*default vx,vy and vz: either observation or 0 */ - if (!iomodel->vx){ - if (iomodel->vx_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts; - else for(i=0;i<6;i++)nodeinputs[i]=0; - this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs)); - this->inputs->AddInput(new PentaVertexInput(VxOldEnum,nodeinputs)); - } - if (!iomodel->vy){ - if (iomodel->vy_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts; - else for(i=0;i<6;i++)nodeinputs[i]=0; - this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs)); - this->inputs->AddInput(new PentaVertexInput(VyOldEnum,nodeinputs)); - } - if (!iomodel->vz){ - if (iomodel->vz_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts; - else for(i=0;i<6;i++)nodeinputs[i]=0; - this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs)); - this->inputs->AddInput(new PentaVertexInput(VzOldEnum,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])); - - //elements of type 3 are macayeal type penta. we collapse the formulation on their base. - if(iomodel->elements_type){ - if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum){ - collapse[analysis_counter]=true; - } - else{ - collapse[analysis_counter]=false; - } - } -} -/*}}}*/ /*Object virtual functions definitions: */ -/*FUNCTION Penta::Echo{{{1*/ - -void Penta::Echo(void){ - this->DeepEcho(); +/*FUNCTION Penta::copy {{{1*/ +Object* Penta::copy() { + + int i; + + Penta* penta=NULL; + + penta=new Penta(); + + /*copy fields: */ + penta->id=this->id; + if(this->inputs){ + penta->inputs=(Inputs*)this->inputs->Copy(); + } + else{ + penta->inputs=new Inputs(); + } + if(this->results){ + penta->results=(Results*)this->results->Copy(); + } + else{ + penta->results=new Results(); + } + /*point parameters: */ + penta->parameters=this->parameters; + + /*now deal with hooks and objects: */ + penta->InitHookNodes(this->numanalyses); + for(i=0;inumanalyses;i++)penta->hnodes[i].copy(&this->hnodes[i]); + penta->hmatice.copy(&this->hmatice); + penta->hmatpar.copy(&this->hmatpar); + penta->hneighbors.copy(&this->hneighbors); + + /*recover objects: */ + penta->nodes=(Node**)xmalloc(6*sizeof(Node*)); //we cannot rely on an analysis_counter to tell us which analysis_type we are running, so we just copy the nodes. + for(i=0;i<6;i++)penta->nodes[i]=this->nodes[i]; + penta->matice=(Matice*)penta->hmatice.delivers(); + penta->matpar=(Matpar*)penta->hmatpar.delivers(); + penta->neighbors=(Penta**)penta->hneighbors.deliverp(); + + return penta; } /*}}}*/ @@ -264,13 +160,66 @@ } /*}}}*/ +/*FUNCTION Penta::Demarshall {{{1*/ +void Penta::Demarshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int i; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*this time, no need to get enum type, the pointer directly points to the beginning of the + *object data (thanks to DataSet::Demarshall):*/ + + memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id); + memcpy(&numanalyses,marshalled_dataset,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses); + + /*allocate dynamic memory: */ + collapse=(bool*)xmalloc(numanalyses*sizeof(bool)); + InitHookNodes(numanalyses); + + /*demarshall hooks: */ + for(i=0;iparameters=NULL; + + /*return: */ + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ +/*FUNCTION Penta::Echo{{{1*/ + +void Penta::Echo(void){ + this->DeepEcho(); +} +/*}}}*/ +/*FUNCTION Penta::Enum {{{1*/ +int Penta::Enum(void){ + + return PentaEnum; + +} +/*}}}*/ /*FUNCTION Penta::Id {{{1*/ int Penta::Id(void){ return id; -} -/*}}}*/ -/*FUNCTION Penta::MyRank {{{1*/ -int Penta::MyRank(void){ - extern int my_rank; - return my_rank; } /*}}}*/ @@ -347,103 +296,288 @@ } /*}}}*/ -/*FUNCTION Penta::Demarshall {{{1*/ -void Penta::Demarshall(char** pmarshalled_dataset){ - - char* marshalled_dataset=NULL; - int i; - - /*recover marshalled_dataset: */ - marshalled_dataset=*pmarshalled_dataset; - - /*this time, no need to get enum type, the pointer directly points to the beginning of the - *object data (thanks to DataSet::Demarshall):*/ - - memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id); - memcpy(&numanalyses,marshalled_dataset,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses); - - /*allocate dynamic memory: */ - collapse=(bool*)xmalloc(numanalyses*sizeof(bool)); - InitHookNodes(numanalyses); - - /*demarshall hooks: */ - for(i=0;iFindParam(&analysis_type,AnalysisTypeEnum); + + /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */ + if (analysis_type==ControlAnalysisEnum){ + InputUpdateFromSolutionDiagnosticHoriz( solution); + } + else if (analysis_type==DiagnosticHorizAnalysisEnum){ + InputUpdateFromSolutionDiagnosticHoriz( solution); + } + else if (analysis_type==DiagnosticStokesAnalysisEnum){ + InputUpdateFromSolutionDiagnosticStokes( solution); + } + else if (analysis_type==SlopeAnalysisEnum){ + InputUpdateFromSolutionSlopeCompute( solution); + } + else if (analysis_type==PrognosticAnalysisEnum){ + InputUpdateFromSolutionPrognostic( solution); + } + else if (analysis_type==Prognostic2AnalysisEnum){ + InputUpdateFromSolutionPrognostic2(solution); + } + else if (analysis_type==BalancedthicknessAnalysisEnum){ + InputUpdateFromSolutionBalancedthickness( solution); + } + else if (analysis_type==Balancedthickness2AnalysisEnum){ + InputUpdateFromSolutionBalancedthickness2( solution); + } + else if (analysis_type==BalancedvelocitiesAnalysisEnum){ + InputUpdateFromSolutionBalancedvelocities( solution); + } + else{ + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type)); + } +} +/*}}}*/ +/*FUNCTION Penta::InputUpdateFromVector(double* vector, int name, int type);{{{1*/ +void Penta::InputUpdateFromVector(double* vector, int name, int type){ + + /*Check that name is an element input*/ + if (!IsInput(name)) return; + + switch(type){ + + case VertexEnum: + + /*New PentaVertexInpu*/ + double values[6]; + + /*Get values on the 6 vertices*/ + for (int i=0;i<6;i++){ + values[i]=vector[this->nodes[i]->GetVertexDof()]; + } + + /*update input*/ + this->inputs->AddInput(new PentaVertexInput(name,values)); + return; + + default: + + ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type)); + } +} +/*}}}*/ +/*FUNCTION Penta::InputUpdateFromVector(int* vector, int name, int type);{{{1*/ +void Penta::InputUpdateFromVector(int* vector, int name, int type){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Penta::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/ +void Penta::InputUpdateFromVector(bool* vector, int name, int type){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ + +/*Element virtual functions definitions: */ +/*FUNCTION Penta::ComputeBasalStress {{{1*/ +void Penta::ComputeBasalStress(Vec sigma_b){ + + int i,j; + const int numgrids=6; + int doflist[numgrids]; + double xyz_list[numgrids][3]; + double xyz_list_tria[3][3]; + + /*Parameters*/ + double rho_ice,gravity; + double surface_normal[3]; + double bed_normal[3]; + double bed; + double basalforce[3]; + double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ + double stresstensor[6]={0.0}; + double viscosity; + int analysis_type,sub_analysis_type; + + int dofv[3]={0,1,2}; + int dofp[1]={3}; + double Jdet2d; + Tria* tria=NULL; + + /*Gauss*/ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + double gauss_weight; + double gauss_coord[4]; + + /*Output*/ + double pressure; + double sigma_xx,sigma_yy,sigma_zz; + double sigma_xy,sigma_xz,sigma_yz; + double surface=0; + double value=0; - /*demarshall inputs and results: */ - inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); - results=(Results*)DataSetDemarshallRaw(&marshalled_dataset); - - /*demarshall internal parameters: */ - memcpy(collapse,marshalled_dataset,numanalyses*sizeof(bool));marshalled_dataset+=numanalyses*sizeof(bool); - - /*parameters: may not exist even yet, so let Configure handle it: */ - this->parameters=NULL; - - /*return: */ - *pmarshalled_dataset=marshalled_dataset; - return; -} -/*}}}*/ -/*FUNCTION Penta::Enum {{{1*/ -int Penta::Enum(void){ - - return PentaEnum; - -} -/*}}}*/ -/*FUNCTION Penta::copy {{{1*/ -Object* Penta::copy() { + /*flags: */ + bool onbed; + + /*parameters: */ + double stokesreconditioning; + + + /*retrive parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum); + + /*Check analysis_types*/ + if (analysis_type!=DiagnosticAnalysisEnum || sub_analysis_type!=StokesAnalysisEnum) ISSMERROR("Not supported yet!"); + + /*recover some inputs: */ + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); + + if(!onbed){ + //put zero + VecSetValue(sigma_b,id-1,0.0,INSERT_VALUES); + return; + } + + /*recovre material parameters: */ + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + for(i=0;i<3;i++){ + for(j=0;j<3;j++){ + xyz_list_tria[i][j]=xyz_list[i][j]; + } + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria(&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights,2); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igGetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + inputs->GetParameterValue(&pressure, &gauss_coord[0],PressureEnum); + + /*Compute Stress*/ + sigma_xx=viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure + sigma_yy=viscosity*epsilon[1]-pressure*stokesreconditioning; + sigma_zz=viscosity*epsilon[2]-pressure*stokesreconditioning; + sigma_xy=viscosity*epsilon[3]; + sigma_xz=viscosity*epsilon[4]; + sigma_yz=viscosity*epsilon[5]; + + /*Get normal vector to the bed */ + SurfaceNormal(&surface_normal[0],xyz_list_tria); + bed_normal[0] = - surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result + bed_normal[1] = - surface_normal[1]; + bed_normal[2] = - surface_normal[2]; + + /*basalforce*/ + basalforce[0] += sigma_xx*bed_normal[0] + sigma_xy*bed_normal[1] + sigma_xz*bed_normal[2]; + basalforce[1] += sigma_xy*bed_normal[0] + sigma_yy*bed_normal[1] + sigma_yz*bed_normal[2]; + basalforce[2] += sigma_xz*bed_normal[0] + sigma_yz*bed_normal[1] + sigma_zz*bed_normal[2]; + + /*Get the Jacobian determinant */ + tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0],gauss_coord); + value+=sigma_zz*Jdet2d*gauss_weight; + surface+=Jdet2d*gauss_weight; + } + value=value/surface; + + /*Add value to output*/ + VecSetValue(sigma_b,id-1,(const double)value,INSERT_VALUES); +} +/*}}}*/ +/*FUNCTION Penta::ComputePressure {{{1*/ +void Penta::ComputePressure(Vec pg){ int i; - - Penta* penta=NULL; - - penta=new Penta(); - - /*copy fields: */ - penta->id=this->id; - if(this->inputs){ - penta->inputs=(Inputs*)this->inputs->Copy(); - } - else{ - penta->inputs=new Inputs(); - } - if(this->results){ - penta->results=(Results*)this->results->Copy(); - } - else{ - penta->results=new Results(); - } - /*point parameters: */ - penta->parameters=this->parameters; - - /*now deal with hooks and objects: */ - penta->InitHookNodes(this->numanalyses); - for(i=0;inumanalyses;i++)penta->hnodes[i].copy(&this->hnodes[i]); - penta->hmatice.copy(&this->hmatice); - penta->hmatpar.copy(&this->hmatpar); - penta->hneighbors.copy(&this->hneighbors); - - /*recover objects: */ - penta->nodes=(Node**)xmalloc(6*sizeof(Node*)); //we cannot rely on an analysis_counter to tell us which analysis_type we are running, so we just copy the nodes. - for(i=0;i<6;i++)penta->nodes[i]=this->nodes[i]; - penta->matice=(Matice*)penta->hmatice.delivers(); - penta->matpar=(Matpar*)penta->hmatpar.delivers(); - penta->neighbors=(Penta**)penta->hneighbors.deliverp(); - - return penta; -} -/*}}}*/ - - -/*Penta functionality: */ + const int numgrids=6; + int doflist[numgrids]; + double pressure[numgrids]; + double rho_ice,g; + double surface[numgrids]; + double xyz_list[numgrids][3]; + double gauss[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + + /*inputs: */ + bool onwater; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + + /*If on water, skip: */ + if(onwater)return; + + /*Get node data: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + + /*pressure is lithostatic: */ + //md.pressure=md.rho_ice*md.g*(md.surface-md.z); a la matlab + + /*Get dof list on which we will plug the pressure values: */ + GetDofList1(&doflist[0]); + + /*recover value of surface at grids: */ + inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum); + + /*pressure is lithostatic: */ + rho_ice=matpar->GetRhoIce(); + g=matpar->GetG(); + for(i=0;iGetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + + /*If on water, return 0: */ + if(onwater)return 0; + + /*Bail out if this element if: + * -> Non collapsed and not on the surface + * -> collapsed (2d model) and not on bed) */ + if ((!collapse && !onsurface) || (collapse && !onbed)){ + return 0; + } + else if (collapse){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + * and compute CostFunction*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + J=tria->CostFunction(); + delete tria; + return J; + } + else{ + + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + J=tria->CostFunction(); + delete tria; + return J; + } +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrix {{{1*/ +void Penta::CreateKMatrix(Mat Kgg){ + + int analysis_type; + + /*retrive parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*if debugging mode, check that all pointers exist*/ + ISSMASSERT(this->nodes && this->matice && this->matpar && this->neighbors && this->parameters && this->inputs); + + /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ + if (analysis_type==ControlAnalysisEnum){ + CreateKMatrixDiagnosticHoriz( Kgg); + } + else if (analysis_type==DiagnosticHorizAnalysisEnum){ + CreateKMatrixDiagnosticHoriz( Kgg); + } + else if (analysis_type==DiagnosticHutterAnalysisEnum){ + CreateKMatrixDiagnosticHutter( Kgg); + } + else if (analysis_type==DiagnosticVertAnalysisEnum){ + CreateKMatrixDiagnosticVert( Kgg); + } + else if (analysis_type==DiagnosticStokesAnalysisEnum){ + CreateKMatrixDiagnosticStokes( Kgg); + } + else if (analysis_type==SlopeAnalysisEnum){ + CreateKMatrixSlopeCompute( Kgg); + } + else if (analysis_type==PrognosticAnalysisEnum){ + CreateKMatrixPrognostic( Kgg); + } + else if (analysis_type==BalancedthicknessAnalysisEnum){ + CreateKMatrixBalancedthickness( Kgg); + } + else if (analysis_type==BalancedvelocitiesAnalysisEnum){ + CreateKMatrixBalancedvelocities( Kgg); + } + else if (analysis_type==ThermalAnalysisEnum){ + CreateKMatrixThermal( Kgg); + } + else if (analysis_type==MeltingAnalysisEnum){ + CreateKMatrixMelting( Kgg); + } + else ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type)); + +} +/*}}}*/ +/*FUNCTION Penta::CreatePVector {{{1*/ +void Penta::CreatePVector(Vec pg){ + + int analysis_type,sub_analysis_type; + + /*retrive parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum); + + /*if debugging mode, check that all pointers exist*/ + ISSMASSERT(this->nodes && this->matice && this->matpar && this->neighbors && this->parameters && this->inputs); + + /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ + if (analysis_type==ControlAnalysisEnum){ + CreatePVectorDiagnosticHoriz( pg); + } + else if (analysis_type==DiagnosticHorizAnalysisEnum){ + CreatePVectorDiagnosticHoriz( pg); + } + else if (analysis_type==DiagnosticHutterAnalysisEnum){ + CreatePVectorDiagnosticHutter( pg); + } + else if (analysis_type==DiagnosticVertAnalysisEnum){ + CreatePVectorDiagnosticVert( pg); + } + else if (analysis_type==DiagnosticStokesAnalysisEnum){ + CreatePVectorDiagnosticStokes( pg); + } + else if (analysis_type==SlopeAnalysisEnum){ + CreatePVectorSlopeCompute( pg); + } + else if (analysis_type==PrognosticAnalysisEnum){ + CreatePVectorPrognostic( pg); + } + else if (analysis_type==BalancedthicknessAnalysisEnum){ + CreatePVectorBalancedthickness( pg); + } + else if (analysis_type==BalancedvelocitiesAnalysisEnum){ + CreatePVectorBalancedvelocities( pg); + } + else if (analysis_type==ThermalAnalysisEnum){ + CreatePVectorThermal( pg); + } + else if (analysis_type==MeltingAnalysisEnum){ + CreatePVectorMelting( pg); + } + else ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type)); + +} +/*}}}*/ +/*FUNCTION Penta::Du {{{1*/ +void Penta::Du(Vec du_g){ + + int i; + Tria* tria=NULL; + + /*inputs: */ + bool onwater; + bool collapse; + bool onsurface; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + + /*If on water, skip: */ + if(onwater)return; + + /*Bail out if this element if: + * -> Non collapsed and not on the surface + * -> collapsed (2d model) and not on bed) */ + if ((!collapse && !onsurface) || (collapse && !onbed)){ + return; + } + else if (collapse){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + * and compute Du*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + tria->Du(du_g); + delete tria; + return; + } + else{ + + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + tria->Du(du_g); + delete tria; + return; + } +} +/*}}}*/ +/*FUNCTION Penta::GetBedList {{{1*/ +void Penta::GetBedList(double* bedlist){ + + const int numgrids=6; + double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + + inputs->GetParameterValues(bedlist,&gaussgrids[0][0],6,BedEnum); + +} +/*}}}*/ +/*FUNCTION Penta::GetMatPar {{{1*/ +void* Penta::GetMatPar(){ + + return matpar; +} +/*}}}*/ +/*FUNCTION Penta::GetNodes {{{1*/ +void Penta::GetNodes(void** vpnodes){ + + int i; + Node** pnodes=NULL; + + /*virtual object: */ + pnodes=(Node**)vpnodes; + + for(i=0;i<6;i++){ + pnodes[i]=nodes[i]; + } +} +/*}}}*/ +/*FUNCTION Penta::GetOnBed {{{1*/ +bool Penta::GetOnBed(){ + + bool onbed; + + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + return onbed; +} +/*}}}*/ +/*FUNCTION Penta::GetShelf {{{1*/ +bool Penta::GetShelf(){ + bool onshelf; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum); + + return onshelf; +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputs(Vec solution){{{1*/ +void Penta::GetSolutionFromInputs(Vec solution){ + + int analysis_type,sub_analysis_type; + + /*retrive parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum); + + + /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */ + if (analysis_type==DiagnosticAnalysisEnum){ + if (sub_analysis_type==HorizAnalysisEnum){ + GetSolutionFromInputsDiagnosticHoriz(solution); + } + else if(sub_analysis_type==VertAnalysisEnum){ + GetSolutionFromInputsDiagnosticVert(solution); + } + else if(sub_analysis_type==StokesAnalysisEnum){ + GetSolutionFromInputsDiagnosticStokes(solution); + } + else ISSMERROR("sub_analysis: %i (%s) not supported yet",sub_analysis_type,EnumAsString(sub_analysis_type)); + } + else{ + ISSMERROR("analysis: %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type)); + } +} +/*}}}*/ +/*FUNCTION Penta::GetThicknessList {{{1*/ +void Penta::GetThicknessList(double* thickness_list){ + + const int numgrids=6; + double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + inputs->GetParameterValues(thickness_list,&gaussgrids[0][0],6,ThicknessEnum); + +} +/*}}}*/ +/*FUNCTION Penta::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/ +void Penta::GetVectorFromInputs(Vec vector,int NameEnum){ + + int i; + const int numvertices=6; + int doflist1[numvertices]; + + /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */ + for(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 Penta::Gradj {{{1*/ +void Penta::Gradj(Vec gradient,int control_type){ + + /*inputs: */ + bool onwater; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + + /*If on water, skip grad (=0): */ + if(onwater)return; + + if (control_type==DragCoefficientEnum){ + GradjDrag(gradient); + } + else if (control_type=RheologyBEnum){ + GradjB(gradient); + } + else ISSMERROR("%s%i","control type not supported yet: ",control_type); +} +/*}}}*/ +/*FUNCTION Penta::GradjDrag {{{1*/ +void Penta::GradjDrag(Vec gradient){ + + int i; + Tria* tria=NULL; + TriaVertexInput* triavertexinput=NULL; + double temp_gradient[6]={0,0,0,0,0,0}; + + /*inputs: */ + bool onwater; + bool onbed; + bool shelf; + int analysis_type,sub_analysis_type; + + /*retrieve parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum); + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + + /*If on water, skip: */ + if(onwater)return; + + /*If on shelf, skip: */ + if(shelf)return; + + /*Bail out if this element does not touch the bedrock: */ + if (!onbed) return; + + if (sub_analysis_type==HorizAnalysisEnum){ + + /*MacAyeal or Pattyn*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->GradjDrag(gradient); + delete tria; + } + else if (sub_analysis_type==StokesAnalysisEnum){ + + /*Stokes*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->GradjDragStokes(gradient); + delete tria; + } + else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); + + +} +/*}}}*/ +/*FUNCTION Penta::InputAXPY(int YEnum, double scalar, int XEnum);{{{1*/ +void Penta::InputAXPY(int YEnum, double scalar, int XEnum){ + + Input* xinput=NULL; + Input* yinput=NULL; + + /*Find x and y inputs: */ + xinput=(Input*)this->inputs->GetInput(XEnum); + yinput=(Input*)this->inputs->GetInput(YEnum); + + /*some checks: */ + if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!"); + if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!"); + + /*Scale: */ + yinput->AXPY(xinput,scalar); +} +/*}}}*/ +/*FUNCTION Penta::InputControlConstrain(int control_type, double cm_min, double cm_max){{{1*/ +void Penta::InputControlConstrain(int control_type, double cm_min, double cm_max){ + + Input* input=NULL; + + /*Find input: */ + input=(Input*)this->inputs->GetInput(control_type); + + /*Do nothing if we don't find it: */ + if(!input)return; + + /*Constrain input using cm_min and cm_max: */ + input->Constrain(cm_min,cm_max); + +} +/*}}}*/ +/*FUNCTION Penta::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/ +void Penta::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){ + + int i; + Input** new_inputs=NULL; + Input** old_inputs=NULL; + int converged=1; + + 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 ",EnumAsString(enums[2*i+0])); + if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(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=0; + } + + /*Assign output pointers:*/ + *pconverged=converged; + +} +/*}}}*/ +/*FUNCTION Penta::InputDepthAverageAtBase{{{1*/ +void Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){ + ISSMERROR("Not implemented yet (see Penta::InputExtrude and Node::FieldDepthAverageAtBase)"); +} +/*}}}*/ +/*FUNCTION Penta::InputDuplicate(int original_enum,int new_enum){{{1*/ +void Penta::InputDuplicate(int original_enum,int new_enum){ + + Input* original=NULL; + Input* copy=NULL; + + /*Make a copy of the original input: */ + original=(Input*)this->inputs->GetInput(original_enum); + copy=(Input*)original->copy(); + + /*Change copy enum to reinitialized_enum: */ + copy->ChangeEnum(new_enum); + + /*Add copy into inputs, it will wipe off the one already there: */ + inputs->AddObject((Input*)copy); +} +/*}}}*/ +/*FUNCTION Penta::InputScale(int enum_type,double scale_factor){{{1*/ +void Penta::InputScale(int enum_type,double scale_factor){ + + Input* input=NULL; + + /*Make a copy of the original input: */ + input=(Input*)this->inputs->GetInput(enum_type); + + /*Scale: */ + input->Scale(scale_factor); +} +/*}}}*/ +/*FUNCTION Penta::InputToResult(int enum_type,int step,double time){{{1*/ +void Penta::InputToResult(int enum_type,int step,double time){ + + int i; + bool found = false; + Input *input = NULL; + + /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ + for (i=0;iinputs->Size();i++){ + input=(Input*)this->inputs->GetObjectByOffset(i); + if (input->EnumType()==enum_type){ + found=true; + break; + } + } + + /*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)); + +} +/*}}}*/ +/*FUNCTION Penta::MassFlux {{{1*/ +double Penta::MassFlux( double* segment){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Penta::MaxAbsVx(double* pmaxabsvx, bool process_units);{{{1*/ +void Penta::MaxAbsVx(double* pmaxabsvx, bool process_units){ + + int i; + int dim; + const int numgrids=6; + double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + double vx_values[numgrids]; + double maxabsvx; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); + + /*now, compute maximum:*/ + maxabsvx=fabs(vx_values[0]); + for(i=1;imaxabsvx)maxabsvx=fabs(vx_values[i]); + } + + /*Assign output pointers:*/ + *pmaxabsvx=maxabsvx; +} +/*}}}*/ +/*FUNCTION Penta::MaxAbsVy(double* pmaxabsvy, bool process_units);{{{1*/ +void Penta::MaxAbsVy(double* pmaxabsvy, bool process_units){ + + int i; + int dim; + const int numgrids=6; + double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + double vy_values[numgrids]; + double maxabsvy; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); + + /*now, compute maximum:*/ + maxabsvy=fabs(vy_values[0]); + for(i=1;imaxabsvy)maxabsvy=fabs(vy_values[i]); + } + + /*Assign output pointers:*/ + *pmaxabsvy=maxabsvy; +} +/*}}}*/ +/*FUNCTION Penta::MaxAbsVz(double* pmaxabsvz, bool process_units);{{{1*/ +void Penta::MaxAbsVz(double* pmaxabsvz, bool process_units){ + + int i; + int dim; + const int numgrids=6; + double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + double vz_values[numgrids]; + double maxabsvz; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); + + /*now, compute maximum:*/ + maxabsvz=fabs(vz_values[0]); + for(i=1;imaxabsvz)maxabsvz=fabs(vz_values[i]); + } + + /*Assign output pointers:*/ + *pmaxabsvz=maxabsvz; +} +/*}}}*/ +/*FUNCTION Penta::MaxVel(double* pmaxvel, bool process_units);{{{1*/ +void Penta::MaxVel(double* pmaxvel, bool process_units){ + + int i; + int dim; + const int numgrids=6; + double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + double vx_values[numgrids]; + double vy_values[numgrids]; + double vz_values[numgrids]; + double vel_values[numgrids]; + double maxvel; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); + inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); + if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); + + /*now, compute maximum of velocity :*/ + if(dim==2){ + for(i=0;imaxvel)maxvel=vel_values[i]; + } + + /*Assign output pointers:*/ + *pmaxvel=maxvel; + +} +/*}}}*/ +/*FUNCTION Penta::MaxVx(double* pmaxvx, bool process_units);{{{1*/ +void Penta::MaxVx(double* pmaxvx, bool process_units){ + + int i; + int dim; + const int numgrids=6; + double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + double vx_values[numgrids]; + double maxvx; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); + + /*now, compute maximum:*/ + maxvx=vx_values[0]; + for(i=1;imaxvx)maxvx=vx_values[i]; + } + + /*Assign output pointers:*/ + *pmaxvx=maxvx; + +} +/*}}}*/ +/*FUNCTION Penta::MaxVy(double* pmaxvy, bool process_units);{{{1*/ +void Penta::MaxVy(double* pmaxvy, bool process_units){ + + int i; + int dim; + const int numgrids=6; + double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + double vy_values[numgrids]; + double maxvy; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); + + /*now, compute maximum:*/ + maxvy=vy_values[0]; + for(i=1;imaxvy)maxvy=vy_values[i]; + } + + /*Assign output pointers:*/ + *pmaxvy=maxvy; + +} +/*}}}*/ +/*FUNCTION Penta::MaxVz(double* pmaxvz, bool process_units);{{{1*/ +void Penta::MaxVz(double* pmaxvz, bool process_units){ + + int i; + int dim; + const int numgrids=6; + double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + double vz_values[numgrids]; + double maxvz; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); + + /*now, compute maximum:*/ + maxvz=vz_values[0]; + for(i=1;imaxvz)maxvz=vz_values[i]; + } + + /*Assign output pointers:*/ + *pmaxvz=maxvz; + +} +/*}}}*/ +/*FUNCTION Penta::MinVel(double* pminvel, bool process_units);{{{1*/ +void Penta::MinVel(double* pminvel, bool process_units){ + + int i; + int dim; + const int numgrids=6; + double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + double vx_values[numgrids]; + double vy_values[numgrids]; + double vz_values[numgrids]; + double vel_values[numgrids]; + double minvel; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); + inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); + if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); + + /*now, compute minimum of velocity :*/ + if(dim==2){ + for(i=0;iFindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); + + /*now, compute minimum:*/ + minvx=vx_values[0]; + for(i=1;iFindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); + + /*now, compute minimum:*/ + minvy=vy_values[0]; + for(i=1;iFindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); + + /*now, compute minimum:*/ + minvz=vz_values[0]; + for(i=1;iGetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + + /*If on water, return 0: */ + if(onwater)return 0; + + /*Bail out if this element if: + * -> Non collapsed and not on the surface + * -> collapsed (2d model) and not on bed) */ + if ((!collapse && !onsurface) || (collapse && !onbed)){ + return 0; + } + else if (collapse){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + * and compute Misfit*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + J=tria->Misfit(); + delete tria; + return J; + } + else{ + + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + J=tria->Misfit(); + delete tria; + return J; + } +} +/*}}}*/ +/*FUNCTION Penta::PatchFill(int* pcount, Patch* patch){{{1*/ +void Penta::PatchFill(int* pcount, Patch* patch){ + + int i; + int count; + int vertices_ids[6]; + + + /*recover pointer: */ + count=*pcount; + + /*will be needed later: */ + for(i=0;i<6;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(count,this->id,vertices_ids,6); + elementresult->PatchFill(count,patch); + + /*increment counter: */ + count++; + } + + /*Assign output pointers:*/ + *pcount=count; +} +/*FUNCTION Penta::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){{{1*/ +void Penta::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){ + + int i; + + /*output: */ + int numrows = 0; + int numvertices = 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? :*/ + numvertices=6; //this is a penta element, with 6 vertices + numnodes=elementresult->NumberOfNodalValues(); //ask result object. + } + + /*Assign output pointers:*/ + *pnumrows=numrows; + *pnumvertices=numvertices; + *pnumnodes=numnodes; + +} +/*}}}*/ +/*FUNCTION Penta::ProcessResultsUnits(void){{{1*/ +void Penta::ProcessResultsUnits(void){ + + int i; + + for(i=0;iresults->Size();i++){ + ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); + elementresult->ProcessUnits(this->parameters); + } + +} +/*}}}*/ +/*FUNCTION Penta::SurfaceArea {{{1*/ +double Penta::SurfaceArea(void){ + + double S; + Tria* tria=NULL; + + /*inputs: */ + bool onwater; + bool collapse; + bool onsurface; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + + /*If on water, return 0: */ + if(onwater)return 0; + + /*Bail out if this element if: + * -> Non collapsed and not on the surface + * -> collapsed (2d model) and not on bed) */ + if ((!collapse && !onsurface) || (collapse && !onbed)){ + return 0; + } + else if (collapse){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + * and compute SurfaceArea*/ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). + S=tria->SurfaceArea(); + delete tria; + return S; + } + else{ + + tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). + S=tria->SurfaceArea(); + delete tria; + return S; + } +} +/*}}}*/ +/*FUNCTION Penta::Update(IoModel* iomodel,int analysis_counter,int analysis_type) {{{1*/ +void Penta::Update(int index,IoModel* iomodel,int analysis_counter,int analysis_type){ + + /*Intermediaries*/ + IssmInt i; + int penta_node_ids[6]; + 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 + } + + /*Recover nodes ids needed to initialize the node hook.*/ + for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook. + penta_node_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab + } + + /*hooks: */ + this->SetHookNodes(penta_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type + + //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]; + 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->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->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->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)); + } + + /*default vx,vy and vz: either observation or 0 */ + if (!iomodel->vx){ + if (iomodel->vx_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts; + else for(i=0;i<6;i++)nodeinputs[i]=0; + this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs)); + this->inputs->AddInput(new PentaVertexInput(VxOldEnum,nodeinputs)); + } + if (!iomodel->vy){ + if (iomodel->vy_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts; + else for(i=0;i<6;i++)nodeinputs[i]=0; + this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs)); + this->inputs->AddInput(new PentaVertexInput(VyOldEnum,nodeinputs)); + } + if (!iomodel->vz){ + if (iomodel->vz_obs) for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts; + else for(i=0;i<6;i++)nodeinputs[i]=0; + this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs)); + this->inputs->AddInput(new PentaVertexInput(VzOldEnum,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])); + + //elements of type 3 are macayeal type penta. we collapse the formulation on their base. + if(iomodel->elements_type){ + if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum){ + collapse[analysis_counter]=true; + } + else{ + collapse[analysis_counter]=false; + } + } +} +/*}}}*/ +/*FUNCTION Penta::UpdateGeometry{{{1*/ +void Penta::UpdateGeometry(void){ + + /*Intermediaries*/ + int i; + double rho_ice,rho_water; + + /*If shelf: hydrostatic equilibrium*/ + if (this->GetShelf()){ + + /*recover material parameters: */ + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoWater(); + + /*Create New Surface: s = (1-rho_ice/rho_water) h*/ + InputDuplicate(ThicknessEnum,SurfaceEnum); //1: copy thickness into surface + InputScale(SurfaceEnum,(1-rho_ice/rho_water)); //2: surface = surface * (1-di) + + /*Create New Bed b = -rho_ice/rho_water h*/ + InputDuplicate(ThicknessEnum,BedEnum); //1: copy thickness into bed + InputScale(BedEnum, -rho_ice/rho_water); //2: bed = bed * (-di) + } + + /*If sheet: surface = bed + thickness*/ + else{ + + /*The bed does not change, update surface only s = b + h*/ + InputDuplicate(BedEnum,SurfaceEnum); //1: copy bed into surface + InputAXPY(SurfaceEnum,1.0,ThicknessEnum); //2: surface = surface + 1 * thickness + } + +} +/*}}}*/ + +/*Penta specific routines: */ /*FUNCTION Penta::SpawnTria {{{1*/ void* Penta::SpawnTria(int g0, int g1, int g2){ @@ -602,40 +1881,4 @@ } /*}}}*/ -/*FUNCTION Penta::InputToResult(int enum_type,int step,double time){{{1*/ -void Penta::InputToResult(int enum_type,int step,double time){ - - int i; - bool found = false; - Input *input = NULL; - - /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */ - for (i=0;iinputs->Size();i++){ - input=(Input*)this->inputs->GetObjectByOffset(i); - if (input->EnumType()==enum_type){ - found=true; - break; - } - } - - /*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)); - -} -/*}}}*/ -/*FUNCTION Penta::ProcessResultsUnits(void){{{1*/ -void Penta::ProcessResultsUnits(void){ - - int i; - - for(i=0;iresults->Size();i++){ - ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); - elementresult->ProcessUnits(this->parameters); - } - -} -/*}}}*/ - -/*updates: */ /*FUNCTION Penta::UpdateFromDakota {{{1*/ void Penta::UpdateFromDakota(void* vinputs){ @@ -645,78 +1888,4 @@ } /*}}}*/ -/*FUNCTION Penta::UpdateGeometry{{{1*/ -void Penta::UpdateGeometry(void){ - - /*Intermediaries*/ - int i; - double rho_ice,rho_water; - - /*If shelf: hydrostatic equilibrium*/ - if (this->GetShelf()){ - - /*recover material parameters: */ - rho_ice=matpar->GetRhoIce(); - rho_water=matpar->GetRhoWater(); - - /*Create New Surface: s = (1-rho_ice/rho_water) h*/ - InputDuplicate(ThicknessEnum,SurfaceEnum); //1: copy thickness into surface - InputScale(SurfaceEnum,(1-rho_ice/rho_water)); //2: surface = surface * (1-di) - - /*Create New Bed b = -rho_ice/rho_water h*/ - InputDuplicate(ThicknessEnum,BedEnum); //1: copy thickness into bed - InputScale(BedEnum, -rho_ice/rho_water); //2: bed = bed * (-di) - } - - /*If sheet: surface = bed + thickness*/ - else{ - - /*The bed does not change, update surface only s = b + h*/ - InputDuplicate(BedEnum,SurfaceEnum); //1: copy bed into surface - InputAXPY(SurfaceEnum,1.0,ThicknessEnum); //2: surface = surface + 1 * thickness - } - -} -/*}}}*/ -/*FUNCTION Penta::InputUpdateFromSolution {{{1*/ -void Penta::InputUpdateFromSolution(double* solution){ - - int analysis_type; - - /*retreive parameters: */ - 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==ControlAnalysisEnum){ - InputUpdateFromSolutionDiagnosticHoriz( solution); - } - else if (analysis_type==DiagnosticHorizAnalysisEnum){ - InputUpdateFromSolutionDiagnosticHoriz( solution); - } - else if (analysis_type==DiagnosticStokesAnalysisEnum){ - InputUpdateFromSolutionDiagnosticStokes( solution); - } - else if (analysis_type==SlopeAnalysisEnum){ - InputUpdateFromSolutionSlopeCompute( solution); - } - else if (analysis_type==PrognosticAnalysisEnum){ - InputUpdateFromSolutionPrognostic( solution); - } - else if (analysis_type==Prognostic2AnalysisEnum){ - InputUpdateFromSolutionPrognostic2(solution); - } - else if (analysis_type==BalancedthicknessAnalysisEnum){ - InputUpdateFromSolutionBalancedthickness( solution); - } - else if (analysis_type==Balancedthickness2AnalysisEnum){ - InputUpdateFromSolutionBalancedthickness2( solution); - } - else if (analysis_type==BalancedvelocitiesAnalysisEnum){ - InputUpdateFromSolutionBalancedvelocities( solution); - } - else{ - ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type)); - } -} -/*Object functions*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHoriz {{{1*/ void Penta::InputUpdateFromSolutionDiagnosticHoriz(double* solution){ @@ -867,32 +2036,4 @@ } /*}}}*/ -/*FUNCTION Penta::GetSolutionFromInputs(Vec solution){{{1*/ -void Penta::GetSolutionFromInputs(Vec solution){ - - int analysis_type,sub_analysis_type; - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum); - - - /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */ - if (analysis_type==DiagnosticAnalysisEnum){ - if (sub_analysis_type==HorizAnalysisEnum){ - GetSolutionFromInputsDiagnosticHoriz(solution); - } - else if(sub_analysis_type==VertAnalysisEnum){ - GetSolutionFromInputsDiagnosticVert(solution); - } - else if(sub_analysis_type==StokesAnalysisEnum){ - GetSolutionFromInputsDiagnosticStokes(solution); - } - else ISSMERROR("sub_analysis: %i (%s) not supported yet",sub_analysis_type,EnumAsString(sub_analysis_type)); - } - else{ - ISSMERROR("analysis: %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type)); - } -} -/*}}}*/ /*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){{{1*/ void Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){ @@ -1002,6 +2143,4 @@ } /*}}}*/ - -/*Object functions*/ /*FUNCTION Penta::IsInput{{{1*/ bool Penta::IsInput(int name){ @@ -1042,278 +2181,4 @@ upper_penta=(Penta*)neighbors[1]; //first one under, second one above return upper_penta; - -} -/*}}}*/ -/*FUNCTION Penta::ComputeBasalStress {{{1*/ -void Penta::ComputeBasalStress(Vec sigma_b){ - - int i,j; - const int numgrids=6; - int doflist[numgrids]; - double xyz_list[numgrids][3]; - double xyz_list_tria[3][3]; - - /*Parameters*/ - double rho_ice,gravity; - double surface_normal[3]; - double bed_normal[3]; - double bed; - double basalforce[3]; - double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ - double stresstensor[6]={0.0}; - double viscosity; - int analysis_type,sub_analysis_type; - - int dofv[3]={0,1,2}; - int dofp[1]={3}; - double Jdet2d; - Tria* tria=NULL; - - /*Gauss*/ - int num_gauss,ig; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* gauss_weights = NULL; - double gauss_weight; - double gauss_coord[4]; - - /*Output*/ - double pressure; - double sigma_xx,sigma_yy,sigma_zz; - double sigma_xy,sigma_xz,sigma_yz; - double surface=0; - double value=0; - - /*flags: */ - bool onbed; - - /*parameters: */ - double stokesreconditioning; - - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum); - - /*Check analysis_types*/ - if (analysis_type!=DiagnosticAnalysisEnum || sub_analysis_type!=StokesAnalysisEnum) ISSMERROR("Not supported yet!"); - - /*recover some inputs: */ - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - - /*retrieve some parameters: */ - this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); - - if(!onbed){ - //put zero - VecSetValue(sigma_b,id-1,0.0,INSERT_VALUES); - return; - } - - /*recovre material parameters: */ - rho_ice=matpar->GetRhoIce(); - gravity=matpar->GetG(); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - for(i=0;i<3;i++){ - for(j=0;j<3;j++){ - xyz_list_tria[i][j]=xyz_list[i][j]; - } - } - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - GaussTria(&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights,2); - - /* Start looping on the number of gaussian points: */ - for (ig=0; igGetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - inputs->GetParameterValue(&pressure, &gauss_coord[0],PressureEnum); - - /*Compute Stress*/ - sigma_xx=viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure - sigma_yy=viscosity*epsilon[1]-pressure*stokesreconditioning; - sigma_zz=viscosity*epsilon[2]-pressure*stokesreconditioning; - sigma_xy=viscosity*epsilon[3]; - sigma_xz=viscosity*epsilon[4]; - sigma_yz=viscosity*epsilon[5]; - - /*Get normal vector to the bed */ - SurfaceNormal(&surface_normal[0],xyz_list_tria); - bed_normal[0] = - surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result - bed_normal[1] = - surface_normal[1]; - bed_normal[2] = - surface_normal[2]; - - /*basalforce*/ - basalforce[0] += sigma_xx*bed_normal[0] + sigma_xy*bed_normal[1] + sigma_xz*bed_normal[2]; - basalforce[1] += sigma_xy*bed_normal[0] + sigma_yy*bed_normal[1] + sigma_yz*bed_normal[2]; - basalforce[2] += sigma_xz*bed_normal[0] + sigma_yz*bed_normal[1] + sigma_zz*bed_normal[2]; - - /*Get the Jacobian determinant */ - tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0],gauss_coord); - value+=sigma_zz*Jdet2d*gauss_weight; - surface+=Jdet2d*gauss_weight; - } - value=value/surface; - - /*Add value to output*/ - VecSetValue(sigma_b,id-1,(const double)value,INSERT_VALUES); -} -/*}}}*/ -/*FUNCTION Penta::ComputePressure {{{1*/ -void Penta::ComputePressure(Vec pg){ - - int i; - const int numgrids=6; - int doflist[numgrids]; - double pressure[numgrids]; - double rho_ice,g; - double surface[numgrids]; - double xyz_list[numgrids][3]; - double gauss[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; - - /*inputs: */ - bool onwater; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*Get node data: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - - /*pressure is lithostatic: */ - //md.pressure=md.rho_ice*md.g*(md.surface-md.z); a la matlab - - /*Get dof list on which we will plug the pressure values: */ - GetDofList1(&doflist[0]); - - /*recover value of surface at grids: */ - inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum); - - /*pressure is lithostatic: */ - rho_ice=matpar->GetRhoIce(); - g=matpar->GetG(); - for(i=0;iGetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&collapse,CollapseEnum); - inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); - - /*If on water, return 0: */ - if(onwater)return 0; - - /*Bail out if this element if: - * -> Non collapsed and not on the surface - * -> collapsed (2d model) and not on bed) */ - if ((!collapse && !onsurface) || (collapse && !onbed)){ - return 0; - } - else if (collapse){ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute CostFunction*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - J=tria->CostFunction(); - delete tria; - return J; - } - else{ - - tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - J=tria->CostFunction(); - delete tria; - return J; - } -} -/*}}}*/ -/*FUNCTION Penta::CreateKMatrix {{{1*/ -void Penta::CreateKMatrix(Mat Kgg){ - - int analysis_type; - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*if debugging mode, check that all pointers exist*/ - ISSMASSERT(this->nodes && this->matice && this->matpar && this->neighbors && this->parameters && this->inputs); - - /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - if (analysis_type==ControlAnalysisEnum){ - CreateKMatrixDiagnosticHoriz( Kgg); - } - else if (analysis_type==DiagnosticHorizAnalysisEnum){ - CreateKMatrixDiagnosticHoriz( Kgg); - } - else if (analysis_type==DiagnosticHutterAnalysisEnum){ - CreateKMatrixDiagnosticHutter( Kgg); - } - else if (analysis_type==DiagnosticVertAnalysisEnum){ - CreateKMatrixDiagnosticVert( Kgg); - } - else if (analysis_type==DiagnosticStokesAnalysisEnum){ - CreateKMatrixDiagnosticStokes( Kgg); - } - else if (analysis_type==SlopeAnalysisEnum){ - CreateKMatrixSlopeCompute( Kgg); - } - else if (analysis_type==PrognosticAnalysisEnum){ - CreateKMatrixPrognostic( Kgg); - } - else if (analysis_type==BalancedthicknessAnalysisEnum){ - CreateKMatrixBalancedthickness( Kgg); - } - else if (analysis_type==BalancedvelocitiesAnalysisEnum){ - CreateKMatrixBalancedvelocities( Kgg); - } - else if (analysis_type==ThermalAnalysisEnum){ - CreateKMatrixThermal( Kgg); - } - else if (analysis_type==MeltingAnalysisEnum){ - CreateKMatrixMelting( Kgg); - } - else ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type)); } @@ -2353,54 +3218,4 @@ } /*}}}*/ -/*FUNCTION Penta::CreatePVector {{{1*/ -void Penta::CreatePVector(Vec pg){ - - int analysis_type,sub_analysis_type; - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum); - - /*if debugging mode, check that all pointers exist*/ - ISSMASSERT(this->nodes && this->matice && this->matpar && this->neighbors && this->parameters && this->inputs); - - /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - if (analysis_type==ControlAnalysisEnum){ - CreatePVectorDiagnosticHoriz( pg); - } - else if (analysis_type==DiagnosticHorizAnalysisEnum){ - CreatePVectorDiagnosticHoriz( pg); - } - else if (analysis_type==DiagnosticHutterAnalysisEnum){ - CreatePVectorDiagnosticHutter( pg); - } - else if (analysis_type==DiagnosticVertAnalysisEnum){ - CreatePVectorDiagnosticVert( pg); - } - else if (analysis_type==DiagnosticStokesAnalysisEnum){ - CreatePVectorDiagnosticStokes( pg); - } - else if (analysis_type==SlopeAnalysisEnum){ - CreatePVectorSlopeCompute( pg); - } - else if (analysis_type==PrognosticAnalysisEnum){ - CreatePVectorPrognostic( pg); - } - else if (analysis_type==BalancedthicknessAnalysisEnum){ - CreatePVectorBalancedthickness( pg); - } - else if (analysis_type==BalancedvelocitiesAnalysisEnum){ - CreatePVectorBalancedvelocities( pg); - } - else if (analysis_type==ThermalAnalysisEnum){ - CreatePVectorThermal( pg); - } - else if (analysis_type==MeltingAnalysisEnum){ - CreatePVectorMelting( pg); - } - else ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type)); - -} -/*}}}*/ /*FUNCTION Penta::CreatePVectorBalancedthickness {{{1*/ void Penta::CreatePVectorBalancedthickness( Vec pg){ @@ -3242,49 +4057,4 @@ } /*}}}*/ -/*FUNCTION Penta::Du {{{1*/ -void Penta::Du(Vec du_g){ - - int i; - Tria* tria=NULL; - - /*inputs: */ - bool onwater; - bool collapse; - bool onsurface; - bool onbed; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&collapse,CollapseEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*Bail out if this element if: - * -> Non collapsed and not on the surface - * -> collapsed (2d model) and not on bed) */ - if ((!collapse && !onsurface) || (collapse && !onbed)){ - return; - } - else if (collapse){ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute Du*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->Du(du_g); - delete tria; - return; - } - else{ - - tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - tria->Du(du_g); - delete tria; - return; - } -} -/*}}}*/ /*FUNCTION Penta::VecExtrude {{{1*/ void Penta::VecExtrude(Vec vector,double* vector_serial,int iscollapsed){ @@ -3389,9 +4159,4 @@ return; -} -/*}}}*/ -/*FUNCTION Penta::InputDepthAverageAtBase{{{1*/ -void Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){ - ISSMERROR("Not implemented yet (see Penta::InputExtrude and Node::FieldDepthAverageAtBase)"); } /*}}}*/ @@ -3558,14 +4323,4 @@ B[i]=dh1dh6[2][i]; } - -} -/*}}}*/ -/*FUNCTION Penta::GetBedList {{{1*/ -void Penta::GetBedList(double* bedlist){ - - const int numgrids=6; - double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; - - inputs->GetParameterValues(bedlist,&gaussgrids[0][0],6,BedEnum); } @@ -4127,10 +4882,4 @@ } /*}}}*/ -/*FUNCTION Penta::GetMatPar {{{1*/ -void* Penta::GetMatPar(){ - - return matpar; -} -/*}}}*/ /*FUNCTION Penta::GetMatrixInvert {{{1*/ void Penta::GetMatrixInvert(double* Ke_invert, double* Ke){ @@ -4376,28 +5125,4 @@ } /*}}}*/ -/*FUNCTION Penta::GetNodes {{{1*/ -void Penta::GetNodes(void** vpnodes){ - - int i; - Node** pnodes=NULL; - - /*virtual object: */ - pnodes=(Node**)vpnodes; - - for(i=0;i<6;i++){ - pnodes[i]=nodes[i]; - } -} -/*}}}*/ -/*FUNCTION Penta::GetOnBed {{{1*/ -bool Penta::GetOnBed(){ - - bool onbed; - - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - - return onbed; -} -/*}}}*/ /*FUNCTION Penta::GetParameterDerivativeValue {{{1*/ void Penta::GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_coord){ @@ -4471,95 +5196,4 @@ } /*}}}*/ -/*FUNCTION Penta::GetShelf {{{1*/ -bool Penta::GetShelf(){ - bool onshelf; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum); - - return onshelf; -} -/*}}}*/ -/*FUNCTION Penta::GetThicknessList {{{1*/ -void Penta::GetThicknessList(double* thickness_list){ - - const int numgrids=6; - double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; - inputs->GetParameterValues(thickness_list,&gaussgrids[0][0],6,ThicknessEnum); - -} -/*}}}*/ -/*FUNCTION Penta::Gradj {{{1*/ -void Penta::Gradj(Vec gradient,int control_type){ - - /*inputs: */ - bool onwater; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - - /*If on water, skip grad (=0): */ - if(onwater)return; - - if (control_type==DragCoefficientEnum){ - GradjDrag(gradient); - } - else if (control_type=RheologyBEnum){ - GradjB(gradient); - } - else ISSMERROR("%s%i","control type not supported yet: ",control_type); -} -/*}}}*/ -/*FUNCTION Penta::GradjDrag {{{1*/ -void Penta::GradjDrag(Vec gradient){ - - int i; - Tria* tria=NULL; - TriaVertexInput* triavertexinput=NULL; - double temp_gradient[6]={0,0,0,0,0,0}; - - /*inputs: */ - bool onwater; - bool onbed; - bool shelf; - int analysis_type,sub_analysis_type; - - /*retrieve parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - parameters->FindParam(&sub_analysis_type,AnalysisTypeEnum); - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*If on shelf, skip: */ - if(shelf)return; - - /*Bail out if this element does not touch the bedrock: */ - if (!onbed) return; - - if (sub_analysis_type==HorizAnalysisEnum){ - - /*MacAyeal or Pattyn*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->GradjDrag(gradient); - delete tria; - } - else if (sub_analysis_type==StokesAnalysisEnum){ - - /*Stokes*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->GradjDragStokes(gradient); - delete tria; - } - else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); - - -} -/*}}}*/ /*FUNCTION Penta::GradjB {{{1*/ void Penta::GradjB(Vec gradient){ @@ -4600,54 +5234,4 @@ -} -/*}}}*/ -/*FUNCTION Penta::MassFlux {{{1*/ -double Penta::MassFlux( double* segment){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Penta::Misfit {{{1*/ -double Penta::Misfit(void){ - - double J; - Tria* tria=NULL; - - /*inputs: */ - bool onwater; - bool collapse; - bool onsurface; - bool onbed; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&collapse,CollapseEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); - - /*If on water, return 0: */ - if(onwater)return 0; - - /*Bail out if this element if: - * -> Non collapsed and not on the surface - * -> collapsed (2d model) and not on bed) */ - if ((!collapse && !onsurface) || (collapse && !onbed)){ - return 0; - } - else if (collapse){ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute Misfit*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - J=tria->Misfit(); - delete tria; - return J; - } - else{ - - tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - J=tria->Misfit(); - delete tria; - return J; - } } /*}}}*/ @@ -4746,49 +5330,4 @@ } /*}}}1*/ -/*FUNCTION Penta::SurfaceArea {{{1*/ -double Penta::SurfaceArea(void){ - - double S; - Tria* tria=NULL; - - /*inputs: */ - bool onwater; - bool collapse; - bool onsurface; - bool onbed; - - /*retrieve inputs :*/ - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&collapse,CollapseEnum); - inputs->GetParameterValue(&onbed,ElementOnBedEnum); - inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); - - /*If on water, return 0: */ - if(onwater)return 0; - - /*Bail out if this element if: - * -> Non collapsed and not on the surface - * -> collapsed (2d model) and not on bed) */ - if ((!collapse && !onsurface) || (collapse && !onbed)){ - return 0; - } - else if (collapse){ - - /*This element should be collapsed into a tria element at its base. Create this tria element, - * and compute SurfaceArea*/ - tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - S=tria->SurfaceArea(); - delete tria; - return S; - } - else{ - - tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - S=tria->SurfaceArea(); - delete tria; - return S; - } -} -/*}}}*/ /*FUNCTION Penta::SurfaceNormal {{{1*/ void Penta::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){ @@ -4817,542 +5356,2 @@ } /*}}}*/ -/*FUNCTION Penta::InputUpdateFromVector(double* vector, int name, int type);{{{1*/ -void Penta::InputUpdateFromVector(double* vector, int name, int type){ - - /*Check that name is an element input*/ - if (!IsInput(name)) return; - - switch(type){ - - case VertexEnum: - - /*New PentaVertexInpu*/ - double values[6]; - - /*Get values on the 6 vertices*/ - for (int i=0;i<6;i++){ - values[i]=vector[this->nodes[i]->GetVertexDof()]; - } - - /*update input*/ - this->inputs->AddInput(new PentaVertexInput(name,values)); - return; - - default: - - ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type)); - } -} -/*}}}*/ -/*FUNCTION Penta::InputUpdateFromVector(int* vector, int name, int type);{{{1*/ -void Penta::InputUpdateFromVector(int* vector, int name, int type){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Penta::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/ -void Penta::InputUpdateFromVector(bool* vector, int name, int type){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Penta::InputUpdateFromConstant(int value, int name);{{{1*/ -void Penta::InputUpdateFromConstant(int constant, int name){ - /*Nothing updated for now*/ -} -/*}}}*/ -/*FUNCTION Penta::InputUpdateFromConstant(bool value, int name);{{{1*/ -void Penta::InputUpdateFromConstant(bool constant, int name){ - /*Nothing updated for now*/ -} -/*}}}*/ -/*FUNCTION Penta::InputUpdateFromConstant(double value, int name);{{{1*/ -void Penta::InputUpdateFromConstant(double constant, int name){ - /*Nothing updated for now*/ -} -/*}}}*/ -/*FUNCTION Penta::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){{{1*/ -void Penta::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){ - - int i; - - /*output: */ - int numrows = 0; - int numvertices = 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? :*/ - numvertices=6; //this is a penta element, with 6 vertices - numnodes=elementresult->NumberOfNodalValues(); //ask result object. - } - - /*Assign output pointers:*/ - *pnumrows=numrows; - *pnumvertices=numvertices; - *pnumnodes=numnodes; - -} -/*}}}*/ -/*FUNCTION Penta::PatchFill(int* pcount, Patch* patch){{{1*/ -void Penta::PatchFill(int* pcount, Patch* patch){ - - int i; - int count; - int vertices_ids[6]; - - - /*recover pointer: */ - count=*pcount; - - /*will be needed later: */ - for(i=0;i<6;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(count,this->id,vertices_ids,6); - elementresult->PatchFill(count,patch); - - /*increment counter: */ - count++; - } - - /*Assign output pointers:*/ - *pcount=count; -} -/*FUNCTION Penta::MinVel(double* pminvel, bool process_units);{{{1*/ -void Penta::MinVel(double* pminvel, bool process_units){ - - int i; - int dim; - const int numgrids=6; - double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; - double vx_values[numgrids]; - double vy_values[numgrids]; - double vz_values[numgrids]; - double vel_values[numgrids]; - double minvel; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); - inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); - if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); - - /*now, compute minimum of velocity :*/ - if(dim==2){ - for(i=0;iFindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); - inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); - if(dim==3) inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); - - /*now, compute maximum of velocity :*/ - if(dim==2){ - for(i=0;imaxvel)maxvel=vel_values[i]; - } - - /*Assign output pointers:*/ - *pmaxvel=maxvel; - -} -/*}}}*/ -/*FUNCTION Penta::MinVx(double* pminvx, bool process_units);{{{1*/ -void Penta::MinVx(double* pminvx, bool process_units){ - - int i; - int dim; - const int numgrids=6; - double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; - double vx_values[numgrids]; - double minvx; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); - - /*now, compute minimum:*/ - minvx=vx_values[0]; - for(i=1;iFindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); - - /*now, compute maximum:*/ - maxvx=vx_values[0]; - for(i=1;imaxvx)maxvx=vx_values[i]; - } - - /*Assign output pointers:*/ - *pmaxvx=maxvx; - -} -/*}}}*/ -/*FUNCTION Penta::MaxAbsVx(double* pmaxabsvx, bool process_units);{{{1*/ -void Penta::MaxAbsVx(double* pmaxabsvx, bool process_units){ - - int i; - int dim; - const int numgrids=6; - double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; - double vx_values[numgrids]; - double maxabsvx; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vx_values[0],&gaussgrids[0][0],numgrids,VxEnum); - - /*now, compute maximum:*/ - maxabsvx=fabs(vx_values[0]); - for(i=1;imaxabsvx)maxabsvx=fabs(vx_values[i]); - } - - /*Assign output pointers:*/ - *pmaxabsvx=maxabsvx; -} -/*}}}*/ -/*FUNCTION Penta::MinVy(double* pminvy, bool process_units);{{{1*/ -void Penta::MinVy(double* pminvy, bool process_units){ - - int i; - int dim; - const int numgrids=6; - double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; - double vy_values[numgrids]; - double minvy; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); - - /*now, compute minimum:*/ - minvy=vy_values[0]; - for(i=1;iFindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); - - /*now, compute maximum:*/ - maxvy=vy_values[0]; - for(i=1;imaxvy)maxvy=vy_values[i]; - } - - /*Assign output pointers:*/ - *pmaxvy=maxvy; - -} -/*}}}*/ -/*FUNCTION Penta::MaxAbsVy(double* pmaxabsvy, bool process_units);{{{1*/ -void Penta::MaxAbsVy(double* pmaxabsvy, bool process_units){ - - int i; - int dim; - const int numgrids=6; - double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; - double vy_values[numgrids]; - double maxabsvy; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vy_values[0],&gaussgrids[0][0],numgrids,VyEnum); - - /*now, compute maximum:*/ - maxabsvy=fabs(vy_values[0]); - for(i=1;imaxabsvy)maxabsvy=fabs(vy_values[i]); - } - - /*Assign output pointers:*/ - *pmaxabsvy=maxabsvy; -} -/*}}}*/ -/*FUNCTION Penta::MinVz(double* pminvz, bool process_units);{{{1*/ -void Penta::MinVz(double* pminvz, bool process_units){ - - int i; - int dim; - const int numgrids=6; - double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; - double vz_values[numgrids]; - double minvz; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); - - /*now, compute minimum:*/ - minvz=vz_values[0]; - for(i=1;iFindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); - - /*now, compute maximum:*/ - maxvz=vz_values[0]; - for(i=1;imaxvz)maxvz=vz_values[i]; - } - - /*Assign output pointers:*/ - *pmaxvz=maxvz; - -} -/*}}}*/ -/*FUNCTION Penta::MaxAbsVz(double* pmaxabsvz, bool process_units);{{{1*/ -void Penta::MaxAbsVz(double* pmaxabsvz, bool process_units){ - - int i; - int dim; - const int numgrids=6; - double gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; - double vz_values[numgrids]; - double maxabsvz; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValues(&vz_values[0],&gaussgrids[0][0],numgrids,VzEnum); - - /*now, compute maximum:*/ - maxabsvz=fabs(vz_values[0]); - for(i=1;imaxabsvz)maxabsvz=fabs(vz_values[i]); - } - - /*Assign output pointers:*/ - *pmaxabsvz=maxabsvz; -} -/*}}}*/ -/*FUNCTION Penta::InputDuplicate(int original_enum,int new_enum){{{1*/ -void Penta::InputDuplicate(int original_enum,int new_enum){ - - Input* original=NULL; - Input* copy=NULL; - - /*Make a copy of the original input: */ - original=(Input*)this->inputs->GetInput(original_enum); - copy=(Input*)original->copy(); - - /*Change copy enum to reinitialized_enum: */ - copy->ChangeEnum(new_enum); - - /*Add copy into inputs, it will wipe off the one already there: */ - inputs->AddObject((Input*)copy); -} -/*}}}*/ -/*FUNCTION Penta::InputScale(int enum_type,double scale_factor){{{1*/ -void Penta::InputScale(int enum_type,double scale_factor){ - - Input* input=NULL; - - /*Make a copy of the original input: */ - input=(Input*)this->inputs->GetInput(enum_type); - - /*Scale: */ - input->Scale(scale_factor); -} -/*}}}*/ -/*FUNCTION Penta::InputAXPY(int YEnum, double scalar, int XEnum);{{{1*/ -void Penta::InputAXPY(int YEnum, double scalar, int XEnum){ - - Input* xinput=NULL; - Input* yinput=NULL; - - /*Find x and y inputs: */ - xinput=(Input*)this->inputs->GetInput(XEnum); - yinput=(Input*)this->inputs->GetInput(YEnum); - - /*some checks: */ - if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!"); - if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!"); - - /*Scale: */ - yinput->AXPY(xinput,scalar); -} -/*}}}*/ -/*FUNCTION Penta::InputControlConstrain(int control_type, double cm_min, double cm_max){{{1*/ -void Penta::InputControlConstrain(int control_type, double cm_min, double cm_max){ - - Input* input=NULL; - - /*Find input: */ - input=(Input*)this->inputs->GetInput(control_type); - - /*Do nothing if we don't find it: */ - if(!input)return; - - /*Constrain input using cm_min and cm_max: */ - input->Constrain(cm_min,cm_max); - -} -/*}}}*/ -/*FUNCTION Penta::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/ -void Penta::GetVectorFromInputs(Vec vector,int NameEnum){ - - int i; - const int numvertices=6; - int doflist1[numvertices]; - - /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */ - for(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 Penta::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/ -void Penta::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){ - - int i; - Input** new_inputs=NULL; - Input** old_inputs=NULL; - int converged=1; - - 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 ",EnumAsString(enums[2*i+0])); - if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(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=0; - } - - /*Assign output pointers:*/ - *pconverged=converged; - -} -/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 4249) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 4250) @@ -62,16 +62,7 @@ void InputUpdateFromConstant(int constant, int name); void InputUpdateFromSolution(double* solutiong); - void InputUpdateFromSolutionBalancedthickness( double* solutiong); - void InputUpdateFromSolutionBalancedthickness2( double* solutiong); - void InputUpdateFromSolutionBalancedvelocities( double* solutiong); - void InputUpdateFromSolutionDiagnosticHoriz( double* solutiong); - void InputUpdateFromSolutionDiagnosticStokes( double* solutiong); - void InputUpdateFromSolutionPrognostic( double* solutiong); - void InputUpdateFromSolutionPrognostic2(double* solutiong); - void InputUpdateFromSolutionSlopeCompute( double* solutiong); void InputUpdateFromVector(bool* vector, int name, int type); void InputUpdateFromVector(double* vector, int name, int type); void InputUpdateFromVector(int* vector, int name, int type); - void UpdateFromDakota(void* inputs); /*}}}*/ /*Element virtual functions definitions: {{{1*/ @@ -178,4 +169,12 @@ Penta* GetUpperElement(void); void InputExtrude(int enum_type); + void InputUpdateFromSolutionBalancedthickness( double* solutiong); + void InputUpdateFromSolutionBalancedthickness2( double* solutiong); + void InputUpdateFromSolutionBalancedvelocities( double* solutiong); + void InputUpdateFromSolutionDiagnosticHoriz( double* solutiong); + void InputUpdateFromSolutionDiagnosticStokes( double* solutiong); + void InputUpdateFromSolutionPrognostic( double* solutiong); + void InputUpdateFromSolutionPrognostic2(double* solutiong); + void InputUpdateFromSolutionSlopeCompute( double* solutiong); bool IsInput(int name); bool IsOnSurface(void); @@ -187,4 +186,5 @@ void* SpawnTria(int g0, int g1, int g2); void SurfaceNormal(double* surface_normal, double xyz_list[3][3]); + void UpdateFromDakota(void* inputs); void VecExtrude(Vec vector,double* vector_serial,int iscollapsed); /*}}}*/