Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 4284) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 4285) @@ -58,4 +58,18 @@ case SlopeAnalysisEnum: + CreateNodesSlopeCompute(pnodes, iomodel,iomodel_handle); + CreateConstraintsSlopeCompute(pconstraints,iomodel,iomodel_handle); + CreateLoadsSlopeCompute(ploads,iomodel,iomodel_handle); + UpdateElementsSlopeCompute(elements,iomodel,iomodel_handle,analysis_counter,analysis_type); + break; + + case BedSlopeAnalysisEnum: + CreateNodesSlopeCompute(pnodes, iomodel,iomodel_handle); + CreateConstraintsSlopeCompute(pconstraints,iomodel,iomodel_handle); + CreateLoadsSlopeCompute(ploads,iomodel,iomodel_handle); + UpdateElementsSlopeCompute(elements,iomodel,iomodel_handle,analysis_counter,analysis_type); + break; + + case SurfaceSlopeAnalysisEnum: CreateNodesSlopeCompute(pnodes, iomodel,iomodel_handle); CreateConstraintsSlopeCompute(pconstraints,iomodel,iomodel_handle); Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4284) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4285) @@ -324,5 +324,5 @@ InputUpdateFromSolutionDiagnosticStokes( solution); } - else if (analysis_type==SlopeAnalysisEnum){ + else if (analysis_type==BedSlopeAnalysisEnum || analysis_type==SurfaceSlopeAnalysisEnum){ InputUpdateFromSolutionSlopeCompute( solution); } @@ -664,5 +664,5 @@ CreateKMatrixDiagnosticStokes( Kgg); } - else if (analysis_type==SlopeAnalysisEnum){ + else if (analysis_type==BedSlopeAnalysisEnum || analysis_type==SurfaceSlopeAnalysisEnum){ CreateKMatrixSlopeCompute( Kgg); } @@ -714,5 +714,5 @@ CreatePVectorDiagnosticStokes( pg); } - else if (analysis_type==SlopeAnalysisEnum){ + else if (analysis_type==BedSlopeAnalysisEnum || analysis_type==SurfaceSlopeAnalysisEnum){ CreatePVectorSlopeCompute( pg); } Index: /issm/trunk/src/c/objects/Elements/Sing.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Sing.cpp (revision 4284) +++ /issm/trunk/src/c/objects/Elements/Sing.cpp (revision 4285) @@ -31,11 +31,54 @@ } /*}}}*/ -/*FUNCTION void Sing::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);{{{1*/ -void Sing::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ /*Object virtual functions definitions: */ +/*FUNCTION Sing::copy {{{1*/ +Object* Sing::copy() { + + int i; + Sing* sing=NULL; + + sing=new Sing(); + + /*copy fields: */ + sing->id=this->id; + if(this->inputs){ + sing->inputs=(Inputs*)this->inputs->Copy(); + } + else{ + sing->inputs=new Inputs(); + } + /*point parameters: */ + sing->parameters=this->parameters; + + /*pointers: */ + sing->node=this->node; + sing->matice=this->matice; + sing->matpar=this->matpar; + + return sing; +} +/*}}}*/ +/*FUNCTION Sing::DeepEcho {{{1*/ +void Sing::DeepEcho(void){ + + printf("Sing:\n"); + printf(" id: %i\n",id); + node->DeepEcho(); + matice->DeepEcho(); + matpar->DeepEcho(); + printf(" parameters\n"); + parameters->DeepEcho(); + printf(" inputs\n"); + inputs->DeepEcho(); + + return; +} +/*}}}*/ +/*FUNCTION Sing::Demarshall {{{1*/ +void Sing::Demarshall(char** pmarshalled_dataset){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ /*FUNCTION Sing::Echo{{{1*/ @@ -53,22 +96,23 @@ } /*}}}*/ -/*FUNCTION Sing::DeepEcho {{{1*/ -void Sing::DeepEcho(void){ - - printf("Sing:\n"); - printf(" id: %i\n",id); - node->DeepEcho(); - matice->DeepEcho(); - matpar->DeepEcho(); - printf(" parameters\n"); - parameters->DeepEcho(); - printf(" inputs\n"); - inputs->DeepEcho(); - - return; +/*FUNCTION Sing::Enum {{{1*/ +int Sing::Enum(void){ + + return SingEnum; + } /*}}}*/ /*FUNCTION Sing::Id {{{1*/ int Sing::Id(void){ return id; } +/*}}}*/ +/*FUNCTION Sing::Marshall {{{1*/ +void Sing::Marshall(char** pmarshalled_dataset){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::MashallSize {{{1*/ +int Sing::MarshallSize(){ + ISSMERROR(" not supported yet!"); +} /*}}}*/ /*FUNCTION Sing::MyRank {{{1*/ @@ -78,55 +122,84 @@ } /*}}}*/ -/*FUNCTION Sing::Marshall {{{1*/ -void Sing::Marshall(char** pmarshalled_dataset){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::MashallSize {{{1*/ -int Sing::MarshallSize(){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::Demarshall {{{1*/ -void Sing::Demarshall(char** pmarshalled_dataset){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::Enum {{{1*/ -int Sing::Enum(void){ - - return SingEnum; - -} -/*}}}*/ -/*FUNCTION Sing::copy {{{1*/ -Object* Sing::copy() { - - int i; - Sing* sing=NULL; - - sing=new Sing(); - - /*copy fields: */ - sing->id=this->id; - if(this->inputs){ - sing->inputs=(Inputs*)this->inputs->Copy(); - } - else{ - sing->inputs=new Inputs(); - } - /*point parameters: */ - sing->parameters=this->parameters; - - /*pointers: */ - sing->node=this->node; - sing->matice=this->matice; - sing->matpar=this->matpar; - - return sing; -} -/*}}}*/ - -/*Sing management*/ + +/*Update virtual functions definitions: */ +/*FUNCTION Sing::InputUpdateFromSolution {{{1*/ +void Sing::InputUpdateFromSolution(double* solution){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::InputUpdateFromVector(double* vector, int name, int type);{{{1*/ +void Sing::InputUpdateFromVector(double* vector, int name, int type){ + + /*Check that name is an element input*/ + if (!IsInput(name)) return; + + switch(type){ + + case VertexEnum: + + /*New SingVertexInpu*/ + double value; + + /*Get values on the 6 vertices*/ + value=vector[node->GetVertexDof()]; + + /*update input*/ + this->inputs->AddInput(new SingVertexInput(name,value)); + return; + + default: + ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type)); + } +} +/*}}}*/ +/*FUNCTION Sing::InputUpdateFromVector(int* vector, int name, int type);{{{1*/ +void Sing::InputUpdateFromVector(int* vector, int name, int type){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/ +void Sing::InputUpdateFromVector(bool* vector, int name, int type){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ + +/*Element virtual functions definitions: */ +/*FUNCTION Sing::ComputeBasalStress {{{1*/ +void Sing::ComputeBasalStress(Vec p_g){ + + ISSMERROR("Not implemented yet"); + +} +/*}}}*/ +/*FUNCTION Sing::ComputePressure {{{1*/ +void Sing::ComputePressure(Vec p_g){ + + int dof; + double pressure; + double thickness; + double rho_ice,g; + + /*Get dof list on which we will plug the pressure values: */ + GetDofList1(&dof); + + /*pressure is lithostatic: */ + rho_ice=matpar->GetRhoIce(); + g=matpar->GetG(); + inputs->GetParameterValue(&thickness,ThicknessEnum); + pressure=rho_ice*g*thickness; + + /*plug local pressure values into global pressure vector: */ + VecSetValue(p_g,dof,pressure,INSERT_VALUES); + +} +/*}}}*/ +/*FUNCTION Sing::ComputeStrainRate {{{1*/ +void Sing::ComputeStrainRate(Vec p_g){ + + ISSMERROR("Not implemented yet"); + +} +/*}}}*/ /*FUNCTION Sing::Configure {{{1*/ void Sing::Configure(Elements* elementsin,Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ @@ -136,4 +209,471 @@ } /*}}}*/ +/*FUNCTION Sing::CostFunction {{{1*/ +double Sing::CostFunction(){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::CreateKMatrix {{{1*/ + +void Sing::CreateKMatrix(Mat Kgg){ + + int analysis_type; + + /*retrive parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ + if (analysis_type==DiagnosticHutterAnalysisEnum){ + CreateKMatrixDiagnosticHutter( Kgg); + } + else{ + ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type)); + } + +} +/*}}}*/ +/*FUNCTION Sing::CreatePVector {{{1*/ +void Sing::CreatePVector(Vec pg){ + + int analysis_type; + + /*retrive parameters: */ + parameters->FindParam(&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 Sing::Du {{{1*/ +void Sing::Du(Vec){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::GetBedList {{{1*/ +void Sing::GetBedList(double*){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::GetMatPar {{{1*/ +void* Sing::GetMatPar(){ + + return matpar; +} +/*}}}*/ +/*FUNCTION Sing::GetNodes {{{1*/ +void Sing::GetNodes(void** vpnodes){ + + Node** pnodes=NULL; + + /*recover nodes: */ + pnodes=(Node**)vpnodes; + + pnodes[0]=node; +} +/*}}}*/ +/*FUNCTION Sing::GetOnBed {{{1*/ +bool Sing::GetOnBed(){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::GetShelf {{{1*/ +bool Sing::GetShelf(){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::GetSolutionFromInputs(Vec solution);{{{1*/ +void Sing::GetSolutionFromInputs(Vec solution){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::GetThicknessList {{{1*/ +void Sing::GetThicknessList(double* thickness_list){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/ +void Sing::GetVectorFromInputs(Vec vector,int NameEnum){ + + int i; + const int numvertices=1; + 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 Sing::Gradj {{{1*/ +void Sing::Gradj(Vec gradient,int control_type){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::GradB {{{1*/ +void Sing::GradjB(Vec gradient){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::GradjDrag {{{1*/ +void Sing::GradjDrag(Vec gradient){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::InputAXPY(int YEnum, double scalar, int XEnum);{{{1*/ +void Sing::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 Sing::InputControlConstrain(int control_type, double cm_min, double cm_max){{{1*/ +void Sing::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 Sing::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/ +void Sing::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 Sing::InputDuplicate(int original_enum,int new_enum){{{1*/ +void Sing::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 Sing::InputScale(int enum_type,double scale_factor){{{1*/ +void Sing::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 Sing::InputToResult(int enum_type,int step,double time){{{1*/ +void Sing::InputToResult(int enum_type,int step,double time){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::MassFlux {{{1*/ +double Sing::MassFlux( double* segment){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::MaxAbsVx(double* pmaxabsvx, bool process_units);{{{1*/ +void Sing::MaxAbsVx(double* pmaxabsvx, bool process_units){ + + int dim; + double maxabsvx; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValue(&maxabsvx,VxEnum); + maxabsvx=fabs(maxabsvx); + + /*Assign output pointers:*/ + *pmaxabsvx=maxabsvx; +} +/*}}}*/ +/*FUNCTION Sing::MaxAbsVy(double* pmaxabsvy, bool process_units);{{{1*/ +void Sing::MaxAbsVy(double* pmaxabsvy, bool process_units){ + + int dim; + double maxabsvy; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValue(&maxabsvy,VyEnum); + maxabsvy=fabs(maxabsvy); + + /*Assign output pointers:*/ + *pmaxabsvy=maxabsvy; +} +/*}}}*/ +/*FUNCTION Sing::MaxAbsVz(double* pmaxabsvz, bool process_units);{{{1*/ +void Sing::MaxAbsVz(double* pmaxabsvz, bool process_units){ + + int dim; + double maxabsvz; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValue(&maxabsvz,VzEnum); + maxabsvz=fabs(maxabsvz); + + /*Assign output pointers:*/ + *pmaxabsvz=maxabsvz; +} +/*}}}*/ +/*FUNCTION Sing::MaxVel(double* pmaxvel, bool process_units);{{{1*/ +void Sing::MaxVel(double* pmaxvel, bool process_units){ + + int dim; + double vx; + double vy; + double vz; + double maxvel; + + /*retrive velocity values at nodes */ + inputs->GetParameterValue(&vx,VxEnum); + inputs->GetParameterValue(&vy,VyEnum); + if(dim==3)inputs->GetParameterValue(&vz,VzEnum); + + /*now, compute maximum of velocity :*/ + if(dim==2){ + maxvel=sqrt(pow(vx,2)+pow(vy,2)); + } + else{ + maxvel=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2)); + } + + /*Assign output pointers:*/ + *pmaxvel=maxvel; + +} +/*}}}*/ +/*FUNCTION Sing::MaxVx(double* pmaxvx, bool process_units);{{{1*/ +void Sing::MaxVx(double* pmaxvx, bool process_units){ + + int dim; + double maxvx; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValue(&maxvx,VxEnum); + + /*Assign output pointers:*/ + *pmaxvx=maxvx; + +} +/*}}}*/ +/*FUNCTION Sing::MaxVy(double* pmaxvy, bool process_units);{{{1*/ +void Sing::MaxVy(double* pmaxvy, bool process_units){ + + int dim; + double maxvy; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValue(&maxvy,VyEnum); + + /*Assign output pointers:*/ + *pmaxvy=maxvy; + +} +/*}}}*/ +/*FUNCTION Sing::MaxVz(double* pmaxvz, bool process_units);{{{1*/ +void Sing::MaxVz(double* pmaxvz, bool process_units){ + + int dim; + double maxvz; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValue(&maxvz,VzEnum); + + /*Assign output pointers:*/ + *pmaxvz=maxvz; + +} +/*}}}*/ +/*FUNCTION Sing::MinVel(double* pminvel, bool process_units);{{{1*/ +void Sing::MinVel(double* pminvel, bool process_units){ + + int dim; + double vx; + double vy; + double vz; + double minvel; + + /*retrive velocity values at nodes */ + inputs->GetParameterValue(&vx,VxEnum); + inputs->GetParameterValue(&vy,VyEnum); + if(dim==3)inputs->GetParameterValue(&vz,VzEnum); + + /*now, compute minimum of velocity :*/ + if(dim==2){ + minvel=sqrt(pow(vx,2)+pow(vy,2)); + } + else{ + minvel=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2)); + } + + /*Assign output pointers:*/ + *pminvel=minvel; + +} +/*}}}*/ +/*FUNCTION Sing::MinVx(double* pminvx, bool process_units);{{{1*/ +void Sing::MinVx(double* pminvx, bool process_units){ + + int dim; + double minvx; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValue(&minvx,VxEnum); + + /*Assign output pointers:*/ + *pminvx=minvx; + +} +/*}}}*/ +/*FUNCTION Sing::MinVy(double* pminvy, bool process_units);{{{1*/ +void Sing::MinVy(double* pminvy, bool process_units){ + + int dim; + double minvy; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValue(&minvy,VyEnum); + + /*Assign output pointers:*/ + *pminvy=minvy; + +} +/*}}}*/ +/*FUNCTION Sing::MinVz(double* pminvz, bool process_units);{{{1*/ +void Sing::MinVz(double* pminvz, bool process_units){ + + int dim; + double minvz; + + /*retrieve dim parameter: */ + parameters->FindParam(&dim,DimEnum); + + /*retrive velocity values at nodes */ + inputs->GetParameterValue(&minvz,VzEnum); + + /*Assign output pointers:*/ + *pminvz=minvz; + +} +/*}}}*/ +/*FUNCTION Sing::Misfit {{{1*/ +double Sing::Misfit(void){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::PatchFill(int* pcount, Patch* patch){{{1*/ +void Sing::PatchFill(int* pcount, Patch* patch){ + + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){{{1*/ +void Sing::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){ + + ISSMERROR(" not supported yet!"); + +} +/*}}}*/ +/*FUNCTION Sing::ProcessResultsUnits(void){{{1*/ +void Sing::ProcessResultsUnits(void){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::SurfaceArea {{{1*/ +double Sing::SurfaceArea( void){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Sing::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);{{{1*/ +void Sing::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ + +/*Sing specific routines: */ /*FUNCTION Sing::IsInput{{{1*/ bool Sing::IsInput(int name){ @@ -143,86 +683,4 @@ } else return false; -} -/*}}}*/ -/*FUNCTION Sing::InputUpdateFromSolution {{{1*/ -void Sing::InputUpdateFromSolution(double* solution){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::GetSolutionFromInputs(Vec solution);{{{1*/ -void Sing::GetSolutionFromInputs(Vec solution){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::InputToResult(int enum_type,int step,double time){{{1*/ -void Sing::InputToResult(int enum_type,int step,double time){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::ProcessResultsUnits(void){{{1*/ -void Sing::ProcessResultsUnits(void){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ - -/*Sing functions*/ -/*FUNCTION Sing::ComputeBasalStress {{{1*/ -void Sing::ComputeBasalStress(Vec p_g){ - - ISSMERROR("Not implemented yet"); - -} -/*}}}*/ -/*FUNCTION Sing::ComputePressure {{{1*/ -void Sing::ComputePressure(Vec p_g){ - - int dof; - double pressure; - double thickness; - double rho_ice,g; - - /*Get dof list on which we will plug the pressure values: */ - GetDofList1(&dof); - - /*pressure is lithostatic: */ - rho_ice=matpar->GetRhoIce(); - g=matpar->GetG(); - inputs->GetParameterValue(&thickness,ThicknessEnum); - pressure=rho_ice*g*thickness; - - /*plug local pressure values into global pressure vector: */ - VecSetValue(p_g,dof,pressure,INSERT_VALUES); - -} -/*}}}*/ -/*FUNCTION Sing::ComputeStrainRate {{{1*/ -void Sing::ComputeStrainRate(Vec p_g){ - - ISSMERROR("Not implemented yet"); - -} -/*}}}*/ -/*FUNCTION Sing::CostFunction {{{1*/ -double Sing::CostFunction(){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::CreateKMatrix {{{1*/ - -void Sing::CreateKMatrix(Mat Kgg){ - - int analysis_type; - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ - if (analysis_type==DiagnosticHutterAnalysisEnum){ - CreateKMatrixDiagnosticHutter( Kgg); - } - else{ - ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type)); - } - } /*}}}*/ @@ -246,22 +704,4 @@ MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES); - -} -/*}}}*/ -/*FUNCTION Sing::CreatePVector {{{1*/ -void Sing::CreatePVector(Vec pg){ - - int analysis_type; - - /*retrive parameters: */ - parameters->FindParam(&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)); - } } @@ -317,14 +757,4 @@ } /*}}}*/ -/*FUNCTION Sing::Du {{{1*/ -void Sing::Du(Vec){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::GetBedList {{{1*/ -void Sing::GetBedList(double*){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ /*FUNCTION Sing::GetDofList {{{1*/ void Sing::GetDofList(int* doflist,int* pnumberofdofspernode){ @@ -352,66 +782,4 @@ } /*}}}*/ -/*FUNCTION Sing::GetMatPar {{{1*/ -void* Sing::GetMatPar(){ - - return matpar; -} -/*}}}*/ -/*FUNCTION Sing::GetNodes {{{1*/ -void Sing::GetNodes(void** vpnodes){ - - Node** pnodes=NULL; - - /*recover nodes: */ - pnodes=(Node**)vpnodes; - - pnodes[0]=node; -} -/*}}}*/ -/*FUNCTION Sing::GetOnBed {{{1*/ -bool Sing::GetOnBed(){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::GetShelf {{{1*/ -bool Sing::GetShelf(){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::GetThicknessList {{{1*/ -void Sing::GetThicknessList(double* thickness_list){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::Gradj {{{1*/ -void Sing::Gradj(Vec gradient,int control_type){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::GradB {{{1*/ -void Sing::GradjB(Vec gradient){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::GradjDrag {{{1*/ -void Sing::GradjDrag(Vec gradient){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::MassFlux {{{1*/ -double Sing::MassFlux( double* segment){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::Misfit {{{1*/ -double Sing::Misfit(void){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::SurfaceArea {{{1*/ -double Sing::SurfaceArea( void){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ /*FUNCTION Sing::SetClone {{{1*/ void Sing::SetClone(int* minranks){ @@ -420,368 +788,2 @@ } /*}}}1*/ -/*FUNCTION Sing::InputUpdateFromVector(double* vector, int name, int type);{{{1*/ -void Sing::InputUpdateFromVector(double* vector, int name, int type){ - - /*Check that name is an element input*/ - if (!IsInput(name)) return; - - switch(type){ - - case VertexEnum: - - /*New SingVertexInpu*/ - double value; - - /*Get values on the 6 vertices*/ - value=vector[node->GetVertexDof()]; - - /*update input*/ - this->inputs->AddInput(new SingVertexInput(name,value)); - return; - - default: - ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type)); - } -} -/*}}}*/ -/*FUNCTION Sing::InputUpdateFromVector(int* vector, int name, int type);{{{1*/ -void Sing::InputUpdateFromVector(int* vector, int name, int type){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/ -void Sing::InputUpdateFromVector(bool* vector, int name, int type){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){{{1*/ -void Sing::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){ - - ISSMERROR(" not supported yet!"); - -} -/*}}}*/ -/*FUNCTION Sing::PatchFill(int* pcount, Patch* patch){{{1*/ -void Sing::PatchFill(int* pcount, Patch* patch){ - - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Sing::MinVel(double* pminvel, bool process_units);{{{1*/ -void Sing::MinVel(double* pminvel, bool process_units){ - - int dim; - double vx; - double vy; - double vz; - double minvel; - - /*retrive velocity values at nodes */ - inputs->GetParameterValue(&vx,VxEnum); - inputs->GetParameterValue(&vy,VyEnum); - if(dim==3)inputs->GetParameterValue(&vz,VzEnum); - - /*now, compute minimum of velocity :*/ - if(dim==2){ - minvel=sqrt(pow(vx,2)+pow(vy,2)); - } - else{ - minvel=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2)); - } - - /*Assign output pointers:*/ - *pminvel=minvel; - -} -/*}}}*/ -/*FUNCTION Sing::MaxVel(double* pmaxvel, bool process_units);{{{1*/ -void Sing::MaxVel(double* pmaxvel, bool process_units){ - - int dim; - double vx; - double vy; - double vz; - double maxvel; - - /*retrive velocity values at nodes */ - inputs->GetParameterValue(&vx,VxEnum); - inputs->GetParameterValue(&vy,VyEnum); - if(dim==3)inputs->GetParameterValue(&vz,VzEnum); - - /*now, compute maximum of velocity :*/ - if(dim==2){ - maxvel=sqrt(pow(vx,2)+pow(vy,2)); - } - else{ - maxvel=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2)); - } - - /*Assign output pointers:*/ - *pmaxvel=maxvel; - -} -/*}}}*/ -/*FUNCTION Sing::MinVx(double* pminvx, bool process_units);{{{1*/ -void Sing::MinVx(double* pminvx, bool process_units){ - - int dim; - double minvx; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValue(&minvx,VxEnum); - - /*Assign output pointers:*/ - *pminvx=minvx; - -} -/*}}}*/ -/*FUNCTION Sing::MaxVx(double* pmaxvx, bool process_units);{{{1*/ -void Sing::MaxVx(double* pmaxvx, bool process_units){ - - int dim; - double maxvx; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValue(&maxvx,VxEnum); - - /*Assign output pointers:*/ - *pmaxvx=maxvx; - -} -/*}}}*/ -/*FUNCTION Sing::MaxAbsVx(double* pmaxabsvx, bool process_units);{{{1*/ -void Sing::MaxAbsVx(double* pmaxabsvx, bool process_units){ - - int dim; - double maxabsvx; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValue(&maxabsvx,VxEnum); - maxabsvx=fabs(maxabsvx); - - /*Assign output pointers:*/ - *pmaxabsvx=maxabsvx; -} -/*}}}*/ -/*FUNCTION Sing::MinVy(double* pminvy, bool process_units);{{{1*/ -void Sing::MinVy(double* pminvy, bool process_units){ - - int dim; - double minvy; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValue(&minvy,VyEnum); - - /*Assign output pointers:*/ - *pminvy=minvy; - -} -/*}}}*/ -/*FUNCTION Sing::MaxVy(double* pmaxvy, bool process_units);{{{1*/ -void Sing::MaxVy(double* pmaxvy, bool process_units){ - - int dim; - double maxvy; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValue(&maxvy,VyEnum); - - /*Assign output pointers:*/ - *pmaxvy=maxvy; - -} -/*}}}*/ -/*FUNCTION Sing::MaxAbsVy(double* pmaxabsvy, bool process_units);{{{1*/ -void Sing::MaxAbsVy(double* pmaxabsvy, bool process_units){ - - int dim; - double maxabsvy; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValue(&maxabsvy,VyEnum); - maxabsvy=fabs(maxabsvy); - - /*Assign output pointers:*/ - *pmaxabsvy=maxabsvy; -} -/*}}}*/ -/*FUNCTION Sing::MinVz(double* pminvz, bool process_units);{{{1*/ -void Sing::MinVz(double* pminvz, bool process_units){ - - int dim; - double minvz; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValue(&minvz,VzEnum); - - /*Assign output pointers:*/ - *pminvz=minvz; - -} -/*}}}*/ -/*FUNCTION Sing::MaxVz(double* pmaxvz, bool process_units);{{{1*/ -void Sing::MaxVz(double* pmaxvz, bool process_units){ - - int dim; - double maxvz; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValue(&maxvz,VzEnum); - - /*Assign output pointers:*/ - *pmaxvz=maxvz; - -} -/*}}}*/ -/*FUNCTION Sing::MaxAbsVz(double* pmaxabsvz, bool process_units);{{{1*/ -void Sing::MaxAbsVz(double* pmaxabsvz, bool process_units){ - - int dim; - double maxabsvz; - - /*retrieve dim parameter: */ - parameters->FindParam(&dim,DimEnum); - - /*retrive velocity values at nodes */ - inputs->GetParameterValue(&maxabsvz,VzEnum); - maxabsvz=fabs(maxabsvz); - - /*Assign output pointers:*/ - *pmaxabsvz=maxabsvz; -} -/*}}}*/ -/*FUNCTION Sing::InputDuplicate(int original_enum,int new_enum){{{1*/ -void Sing::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 Sing::InputScale(int enum_type,double scale_factor){{{1*/ -void Sing::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 Sing::InputAXPY(int YEnum, double scalar, int XEnum);{{{1*/ -void Sing::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 Sing::InputControlConstrain(int control_type, double cm_min, double cm_max){{{1*/ -void Sing::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 Sing::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/ -void Sing::GetVectorFromInputs(Vec vector,int NameEnum){ - - int i; - const int numvertices=1; - 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 Sing::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/ -void Sing::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/Sing.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Sing.h (revision 4284) +++ /issm/trunk/src/c/objects/Elements/Sing.h (revision 4285) @@ -42,22 +42,22 @@ /*}}}*/ /*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 Echo(); + int Id(); + int MyRank(); + void Marshall(char** pmarshalled_dataset); + int MarshallSize(); + int Enum(); + void DeepEcho(); + void Demarshall(char** pmarshalled_dataset); /*}}}*/ /*Update virtual functions resolution: {{{1*/ + void InputUpdateFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");} + void InputUpdateFromConstant(double constant, int name){ISSMERROR("Not implemented yet!");} + void InputUpdateFromConstant(int constant, int name){ISSMERROR("Not implemented yet!");} + void InputUpdateFromSolution(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 InputUpdateFromVector(bool* vector, int name, int type); - void InputUpdateFromConstant(double constant, int name){ISSMERROR("Not implemented yet!");} - void InputUpdateFromConstant(int constant, int name){ISSMERROR("Not implemented yet!");} - void InputUpdateFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");} - void InputUpdateFromSolution(double* solutiong); /*}}}*/ /*Element virtual functions definitions: {{{1*/ Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4284) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4285) @@ -69,161 +69,119 @@ } /*}}}*/ -/*FUNCTION Tria::Update(IoModel* iomodel,int analysis_counter,int analysis_type){{{1*/ -void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index - - /*Intermediaries*/ - int i; - int tria_node_ids[3]; - int tria_vertex_ids[3]; - double nodeinputs[3]; - - /*Checks if debuging*/ - /*{{{2*/ - ISSMASSERT(iomodel->elements); - /*}}}*/ - - /*Recover vertices ids needed to initialize inputs*/ - for(i=0;i<3;i++){ - tria_vertex_ids[i]=(int)iomodel->elements[3*index+i]; //ids for vertices are in the elements array from Matlab - } - - /*Recover nodes ids needed to initialize the node hook.*/ - if (analysis_type==Prognostic2AnalysisEnum || analysis_type==Balancedthickness2AnalysisEnum){ - /*Discontinuous Galerkin*/ - tria_node_ids[0]=iomodel->nodecounter+3*index+1; - tria_node_ids[1]=iomodel->nodecounter+3*index+2; - tria_node_ids[2]=iomodel->nodecounter+3*index+3; + +/*Object virtual functions definitions:*/ +/*FUNCTION Tria::copy {{{1*/ +Object* Tria::copy() { + + int i; + Tria* tria=NULL; + + tria=new Tria(); + + /*copy fields: */ + tria->id=this->id; + tria->interpolation_type=this->interpolation_type; + if(this->inputs){ + tria->inputs=(Inputs*)this->inputs->Copy(); } else{ - /*Continuous Galerkin*/ - for(i=0;i<3;i++){ - tria_node_ids[i]=iomodel->nodecounter+(int)*(iomodel->elements+3*index+i); //ids for vertices are in the elements array from Matlab - } - } - - /*hooks: */ - this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type + tria->inputs=new Inputs(); + } + if(this->results){ + tria->results=(Results*)this->results->Copy(); + } + else{ + tria->results=new Results(); + } + /*point parameters: */ + tria->parameters=this->parameters; + + /*now deal with hooks and objects: */ + tria->InitHookNodes(this->numanalyses); + for(i=0;inumanalyses;i++)tria->hnodes[i].copy(&this->hnodes[i]); + tria->hmatice.copy(&this->hmatice); + tria->hmatpar.copy(&this->hmatpar); + + /*recover objects: */ + tria->nodes=(Node**)xmalloc(3*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<3;i++)tria->nodes[i]=this->nodes[i]; + tria->matice=(Matice*)tria->hmatice.delivers(); + tria->matpar=(Matpar*)tria->hmatpar.delivers(); + + return tria; +} +/*}}}*/ +/*FUNCTION Tria::DeepEcho{{{1*/ +void Tria::DeepEcho(void){ + + printf("Tria:\n"); + printf(" id: %i\n",id); + if(nodes){ + nodes[0]->DeepEcho(); + nodes[1]->DeepEcho(); + nodes[2]->DeepEcho(); + } + else printf("nodes = NULL\n"); + + if (matice) matice->DeepEcho(); + else printf("matice = NULL\n"); + + if (matpar) matpar->DeepEcho(); + else printf("matpar = NULL\n"); + + printf(" parameters\n"); + if (parameters) parameters->DeepEcho(); + else printf("parameters = NULL\n"); + + printf(" inputs\n"); + if (inputs) inputs->DeepEcho(); + else printf("inputs=NULL\n"); + + if (results) results->DeepEcho(); + else printf("results=NULL\n"); - /*add as many inputs per element as requested:*/ - if (iomodel->thickness) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,nodeinputs)); - } - if (iomodel->surface) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,nodeinputs)); - } - if (iomodel->bed) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(BedEnum,nodeinputs)); - } - if (iomodel->drag_coefficient) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->drag_coefficient[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(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<3;i++)nodeinputs[i]=iomodel->melting_rate[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(MeltingRateEnum,nodeinputs)); - } - if (iomodel->accumulation_rate) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(AccumulationRateEnum,nodeinputs)); - } - if (iomodel->geothermalflux) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->geothermalflux[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(GeothermalFluxEnum,nodeinputs)); - } - if (iomodel->dhdt) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->dhdt[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(DhDtEnum,nodeinputs)); - } - if (iomodel->pressure) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->pressure[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(PressureEnum,nodeinputs)); - } - if (iomodel->temperature) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->temperature[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(TemperatureEnum,nodeinputs)); - } - /*vx,vy and vz: */ - if (iomodel->vx) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs)); - this->inputs->AddInput(new TriaVertexInput(VxOldEnum,nodeinputs)); - } - if (iomodel->vy) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs)); - this->inputs->AddInput(new TriaVertexInput(VyOldEnum,nodeinputs)); - } - if (iomodel->vz) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs)); - this->inputs->AddInput(new TriaVertexInput(VzOldEnum,nodeinputs)); - } - if (iomodel->vx_obs) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(VxObsEnum,nodeinputs)); - } - if (iomodel->vy_obs) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(VyObsEnum,nodeinputs)); - } - if (iomodel->vz_obs) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(VzObsEnum,nodeinputs)); - } - if (iomodel->weights) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->weights[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(WeightsEnum,nodeinputs)); - } - if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index])); - if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index])); - if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index])); - if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index])); - - /*Defaults if not provided in iomodel*/ - switch(analysis_type){ - - case DiagnosticHorizAnalysisEnum: case DiagnosticVertAnalysisEnum: case DiagnosticStokesAnalysisEnum: - - /*default vx,vy and vz: either observation or 0 */ - if(!iomodel->vx){ - if (iomodel->vx_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts; - else for(i=0;i<3;i++)nodeinputs[i]=0; - this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs)); - this->inputs->AddInput(new TriaVertexInput(VxOldEnum,nodeinputs)); - } - if(!iomodel->vy){ - if (iomodel->vy_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts; - else for(i=0;i<3;i++)nodeinputs[i]=0; - this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs)); - this->inputs->AddInput(new TriaVertexInput(VyOldEnum,nodeinputs)); - } - if(!iomodel->vz){ - if (iomodel->vz_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_vertex_ids[i]-1]/iomodel->yts; - else for(i=0;i<3;i++)nodeinputs[i]=0; - this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs)); - this->inputs->AddInput(new TriaVertexInput(VzOldEnum,nodeinputs)); - } - break; - - default: - /*No update for other solution types*/ - break; - - } - - //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. + return; +} +/*}}}*/ +/*FUNCTION Tria::Demarshall {{{1*/ +void Tria::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(&interpolation_type,marshalled_dataset,sizeof(interpolation_type));marshalled_dataset+=sizeof(interpolation_type); + memcpy(&numanalyses,marshalled_dataset,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses); + + /*allocate dynamic memory: */ + InitHookNodes(numanalyses); + + /*demarshall hooks: */ + for(i=0;iparameters=NULL; -} -/*}}}*/ - -/*Object virtual functions definitions:*/ + /*return: */ + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ /*FUNCTION Tria::Echo{{{1*/ void Tria::Echo(void){ @@ -255,44 +213,13 @@ } /*}}}*/ -/*FUNCTION Tria::DeepEcho{{{1*/ -void Tria::DeepEcho(void){ - - printf("Tria:\n"); - printf(" id: %i\n",id); - if(nodes){ - nodes[0]->DeepEcho(); - nodes[1]->DeepEcho(); - nodes[2]->DeepEcho(); - } - else printf("nodes = NULL\n"); - - if (matice) matice->DeepEcho(); - else printf("matice = NULL\n"); - - if (matpar) matpar->DeepEcho(); - else printf("matpar = NULL\n"); - - printf(" parameters\n"); - if (parameters) parameters->DeepEcho(); - else printf("parameters = NULL\n"); - - printf(" inputs\n"); - if (inputs) inputs->DeepEcho(); - else printf("inputs=NULL\n"); - - if (results) results->DeepEcho(); - else printf("results=NULL\n"); - - return; +/*FUNCTION Tria::Enum {{{1*/ +int Tria::Enum(void){ + + return TriaEnum; + } /*}}}*/ /*FUNCTION Tria::Id {{{1*/ int Tria::Id(){ return id; } -/*}}}*/ -/*FUNCTION Tria::MyRank {{{1*/ -int Tria::MyRank(void){ - extern int my_rank; - return my_rank; -} /*}}}*/ /*FUNCTION Tria::Marshall {{{1*/ @@ -366,271 +293,25 @@ } /*}}}*/ -/*FUNCTION Tria::Demarshall {{{1*/ -void Tria::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(&interpolation_type,marshalled_dataset,sizeof(interpolation_type));marshalled_dataset+=sizeof(interpolation_type); - memcpy(&numanalyses,marshalled_dataset,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses); - - /*allocate dynamic memory: */ - InitHookNodes(numanalyses); - - /*demarshall hooks: */ - for(i=0;iparameters=NULL; - - /*return: */ - *pmarshalled_dataset=marshalled_dataset; - return; -} -/*}}}*/ -/*FUNCTION Tria::Enum {{{1*/ -int Tria::Enum(void){ - - return TriaEnum; - -} -/*}}}*/ -/*FUNCTION Tria::copy {{{1*/ -Object* Tria::copy() { - - int i; - Tria* tria=NULL; - - tria=new Tria(); - - /*copy fields: */ - tria->id=this->id; - tria->interpolation_type=this->interpolation_type; - if(this->inputs){ - tria->inputs=(Inputs*)this->inputs->Copy(); - } - else{ - tria->inputs=new Inputs(); - } - if(this->results){ - tria->results=(Results*)this->results->Copy(); - } - else{ - tria->results=new Results(); - } - /*point parameters: */ - tria->parameters=this->parameters; - - /*now deal with hooks and objects: */ - tria->InitHookNodes(this->numanalyses); - for(i=0;inumanalyses;i++)tria->hnodes[i].copy(&this->hnodes[i]); - tria->hmatice.copy(&this->hmatice); - tria->hmatpar.copy(&this->hmatpar); - - /*recover objects: */ - tria->nodes=(Node**)xmalloc(3*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<3;i++)tria->nodes[i]=this->nodes[i]; - tria->matice=(Matice*)tria->hmatice.delivers(); - tria->matpar=(Matpar*)tria->hmatpar.delivers(); - - return tria; -} -/*}}}*/ - -/*Tria management: */ -/*FUNCTION Tria::Configure {{{1*/ -void Tria::Configure(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ - - int analysis_counter; - - /*go into parameters and get the analysis_counter: */ - parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); - - /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective - * datasets, using internal ids and offsets hidden in hooks: */ - this->hnodes[analysis_counter].configure(nodesin); - this->hmatice.configure(materialsin); - this->hmatpar.configure(materialsin); - - /*Now, go pick up the objects inside the hooks: */ - this->nodes=(Node**)this->hnodes[analysis_counter].deliverp(); - this->matice=(Matice*)this->hmatice.delivers(); - this->matpar=(Matpar*)this->hmatpar.delivers(); - - /*point parameters to real dataset: */ - this->parameters=parametersin; - -} -/*}}}*/ -/*FUNCTION Tria::SpawnBeam {{{1*/ -void* Tria::SpawnBeam(int g0, int g1){ - - int i; - - /*out of grids g0,g1 and g2 from Tria, build a beam element: */ - Beam* beam=NULL; - int indices[2]; - int zero=0; - Parameters *beam_parameters = NULL; - Inputs *beam_inputs = NULL; - - indices[0]=g0; - indices[1]=g1; - - beam_parameters=this->parameters; - beam_inputs=(Inputs*)this->inputs->SpawnBeamInputs(indices); - - beam=new Beam(); - beam->id=this->id; - beam->inputs=beam_inputs; - beam->parameters=beam_parameters; - - /*now deal with nodes, matice and matpar: */ - beam->nodes=(Node**)xmalloc(2*sizeof(Node*)); - for(i=0;i<2;i++)beam->nodes[i]=this->nodes[indices[i]]; - beam->matice=this->matice; - beam->matpar=this->matpar; - - - return beam; -} -/*}}}*/ -/*FUNCTION Tria::SpawnSing {{{1*/ -void* Tria::SpawnSing(int index){ - - Sing* sing=NULL; - int zero=0; - Parameters *sing_parameters = NULL; - Inputs *sing_inputs = NULL; - - sing_parameters=this->parameters; - sing_inputs=(Inputs*)this->inputs->SpawnSingInputs(index); - - sing=new Sing(); - sing->id=this->id; - sing->inputs=sing_inputs; - sing->parameters=sing_parameters; - - /*now deal with node,matice and matpar: */ - sing->node=this->nodes[index]; - sing->matice=this->matice; - sing->matpar=this->matpar; - - return sing; -} -/*}}}*/ -/*FUNCTION Tria::InputToResult(int enum_type,int step,double time){{{1*/ -void Tria::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 Tria::ProcessResultsUnits(void){{{1*/ -void Tria::ProcessResultsUnits(void){ - - int i; - - for(i=0;iresults->Size();i++){ - ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); - elementresult->ProcessUnits(this->parameters); - } - -} -/*}}}*/ - -/*Updates: */ -/*FUNCTION Tria::UpdateGeometry{{{1*/ -void Tria::UpdateGeometry(void){ - - /*Intermediaries*/ - 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 Tria::UpdateFromDakota {{{1*/ -void Tria::UpdateFromDakota(void* vinputs){ - - int i; - int dofs[1]={0}; - double temperature_list[3]; - double temperature_average; - double B_list[3]; - double B_average; - double new_h[3]; - - /*Update internal data if inputs holds new values: */ - /*inputs->Recover("thickness",&this->properties.h[0],1,dofs,3,(void**)nodes); - if(inputs->Recover("thickness",&new_h[0],1,dofs,3,(void**)nodes)){ - //density, needed later: - double di=(this->matpar->GetRhoIce()/this->matpar->GetRhoWater()); - //Go through grids: - for (i=0;i<3;i++){ - if(nodes[i]->IsOnShelf()){ - this->b[i]=this->b[i]-di*(new_h[i]-h[i]); //hydrostatic equilibrium; - } - this->s[i]=this->b[i]+new_h[i]; - this->h[i]=new_h[i]; - } - }*/ - - ISSMERROR("not supported yet!"); - +/*FUNCTION Tria::MyRank {{{1*/ +int Tria::MyRank(void){ + extern int my_rank; + return my_rank; +} +/*}}}*/ + +/*Update virtual functions definitions: */ +/*FUNCTION Tria::InputUpdateFromConstant(int value, int name);{{{1*/ +void Tria::InputUpdateFromConstant(int constant, int name){ + /*Nothing updated for now*/ +} +/*}}}*/ +/*FUNCTION Tria::InputUpdateFromConstant(double value, int name);{{{1*/ +void Tria::InputUpdateFromConstant(double constant, int name){ + /*Nothing updated for now*/ +} +/*}}}*/ +/*FUNCTION Tria::InputUpdateFromConstant(bool value, int name);{{{1*/ +void Tria::InputUpdateFromConstant(bool constant, int name){ + /*Nothing updated for now*/ } /*}}}*/ @@ -650,5 +331,5 @@ InputUpdateFromSolutionDiagnosticHoriz( solution); } - else if (analysis_type==SlopeAnalysisEnum){ + else if (analysis_type==BedSlopeAnalysisEnum || analysis_type==SurfaceSlopeAnalysisEnum){ InputUpdateFromSolutionSlopeCompute( solution); } @@ -673,220 +354,44 @@ } /*}}}*/ -/*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHoriz {{{1*/ -void Tria::InputUpdateFromSolutionDiagnosticHoriz(double* solution){ - - int i; - - const int numvertices=3; - const int numdofpervertex=2; - const int numdof=numdofpervertex*numvertices; - - int doflist[numdof]; - double values[numdof]; - double vx[numvertices]; - double vy[numvertices]; - double vz[numvertices]; - double vel[numvertices]; - double pressure[numvertices]; - double thickness[numvertices]; - double rho_ice,g; - double gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; - - int dummy; - Input* VzInput=NULL; - double* VzPtr=NULL; - - /*Get dof list: */ - GetDofList(&doflist[0],&dummy); - - /*Use the dof list to index into the solution vector: */ - for(i=0;iGetInput(VzEnum); - if (VzInput){ - if (VzInput->Enum()!=TriaVertexInputEnum){ - ISSMERROR("Cannot compute Vel as Vz is of type %s",EnumAsString(VzInput->Enum())); - } - VzInput->GetValuesPtr(&VzPtr,&dummy); - for(i=0;iGetRhoIce(); - g=matpar->GetG(); - inputs->GetParameterValues(&thickness[0],&gauss[0][0],3,ThicknessEnum); - - for(i=0;iinputs->ChangeEnum(VxEnum,VxOldEnum); - this->inputs->ChangeEnum(VyEnum,VyOldEnum); - this->inputs->ChangeEnum(PressureEnum,PressureOldEnum); - - /*Add vx and vy as inputs to the tria element: */ - this->inputs->AddInput(new TriaVertexInput(VxEnum,vx)); - this->inputs->AddInput(new TriaVertexInput(VyEnum,vy)); - this->inputs->AddInput(new TriaVertexInput(VelEnum,vel)); - this->inputs->AddInput(new TriaVertexInput(PressureEnum,pressure)); - -} - -/*}}}*/ -/*FUNCTION Tria::InputUpdateFromSolutionSlopeCompute {{{1*/ -void Tria::InputUpdateFromSolutionSlopeCompute(double* solution){ +/*FUNCTION Tria::InputUpdateFromVector(double* vector, int name, int type);{{{1*/ +void Tria::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[3]; + + /*Get values on the 6 vertices*/ + for (int i=0;i<3;i++){ + values[i]=vector[this->nodes[i]->GetVertexDof()]; + } + + /*update input*/ + this->inputs->AddInput(new TriaVertexInput(name,values)); + return; + + default: + + ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type)); + } +} +/*}}}*/ +/*FUNCTION Tria::InputUpdateFromVector(int* vector, int name, int type);{{{1*/ +void Tria::InputUpdateFromVector(int* vector, int name, int type){ ISSMERROR(" not supported yet!"); } /*}}}*/ -/*FUNCTION Tria::InputUpdateFromSolutionPrognostic {{{1*/ -void Tria::InputUpdateFromSolutionPrognostic(double* solution){ - - int i; - - const int numvertices=3; - const int numdofpervertex=1; - const int numdof=numdofpervertex*numvertices; - - int doflist[numdof]; - double values[numdof]; - double thickness[numvertices]; - - int dummy; - - /*Get dof list: */ - GetDofList(&doflist[0],&dummy); - - /*Use the dof list to index into the solution vector: */ - for(i=0;iinputs->AddInput(new TriaVertexInput(ThicknessEnum,values)); -} -/*}}}*/ -/*FUNCTION Tria::InputUpdateFromSolutionPrognostic2 {{{1*/ -void Tria::InputUpdateFromSolutionPrognostic2(double* solution){ +/*FUNCTION Tria::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/ +void Tria::InputUpdateFromVector(bool* vector, int name, int type){ ISSMERROR(" not supported yet!"); } /*}}}*/ -/*FUNCTION Tria::InputUpdateFromSolutionBalancedthickness {{{1*/ -void Tria::InputUpdateFromSolutionBalancedthickness(double* solution){ - - int i; - - const int numvertices=3; - const int numdofpervertex=1; - const int numdof=numdofpervertex*numvertices; - - int doflist[numdof]; - double values[numdof]; - double thickness[numvertices]; - - int dummy; - - /*Get dof list: */ - GetDofList(&doflist[0],&dummy); - - /*Use the dof list to index into the solution vector: */ - for(i=0;iinputs->AddInput(new TriaVertexInput(ThicknessEnum,values)); -} -/*}}}*/ -/*FUNCTION Tria::InputUpdateFromSolutionBalancedthickness2 {{{1*/ -void Tria::InputUpdateFromSolutionBalancedthickness2(double* solution){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Tria::InputUpdateFromSolutionBalancedvelocities {{{1*/ -void Tria::InputUpdateFromSolutionBalancedvelocities(double* solution){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Tria::GetSolutionFromInputs(Vec solution){{{1*/ -void Tria::GetSolutionFromInputs(Vec solution){ - - int analysis_type; - - /*retrive 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==DiagnosticHorizAnalysisEnum) - GetSolutionFromInputsDiagnosticHoriz(solution); - else - ISSMERROR("%s%i%s\n","analysis: ",analysis_type," not supported yet"); - -} -/*}}}*/ -/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution){{{1*/ -void Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution){ - - int i; - - const int numvertices=3; - const int numdofpervertex=2; - const int numdof=numdofpervertex*numvertices; - double gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; - - int doflist[numdof]; - double values[numdof]; - double vx; - double vy; - - int dummy; - - /*Get dof list: */ - GetDofList(&doflist[0],&dummy); - - /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ - /*P1 element only for now*/ - for(i=0;iGetParameterValue(&vx,&gauss[i][0],VxEnum); - inputs->GetParameterValue(&vy,&gauss[i][0],VyEnum); - values[i*numdofpervertex+0]=vx; - values[i*numdofpervertex+1]=vy; - } - - /*Add value to global vector*/ - VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); - -} -/*}}}*/ - -/*Object functions*/ -/*FUNCTION Tria::IsInput{{{1*/ -bool Tria::IsInput(int name){ - if (name==SurfaceSlopeXEnum || - name==SurfaceSlopeYEnum){ - return true; - } - else return false; -} -/*}}}*/ + +/*Element virtual functions definitions: */ /*FUNCTION Tria::ComputeBasalStress {{{1*/ void Tria::ComputeBasalStress(Vec eps){ @@ -947,4 +452,28 @@ } /*}}}*/ +/*FUNCTION Tria::Configure {{{1*/ +void Tria::Configure(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ + + int analysis_counter; + + /*go into parameters and get the analysis_counter: */ + parametersin->FindParam(&analysis_counter,AnalysisCounterEnum); + + /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective + * datasets, using internal ids and offsets hidden in hooks: */ + this->hnodes[analysis_counter].configure(nodesin); + this->hmatice.configure(materialsin); + this->hmatpar.configure(materialsin); + + /*Now, go pick up the objects inside the hooks: */ + this->nodes=(Node**)this->hnodes[analysis_counter].deliverp(); + this->matice=(Matice*)this->hmatice.delivers(); + this->matpar=(Matpar*)this->hmatpar.delivers(); + + /*point parameters to real dataset: */ + this->parameters=parametersin; + +} +/*}}}*/ /*FUNCTION Tria::CostFunction {{{1*/ double Tria::CostFunction(void){ @@ -1070,5 +599,5 @@ CreateKMatrixDiagnosticHutter( Kgg); } - else if (analysis_type==SlopeAnalysisEnum){ + else if (analysis_type==BedSlopeAnalysisEnum || analysis_type==SurfaceSlopeAnalysisEnum){ CreateKMatrixSlopeCompute( Kgg); } @@ -1091,1368 +620,4 @@ ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type)); } - -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixBalancedthickness {{{1*/ -void Tria::CreateKMatrixBalancedthickness(Mat Kgg){ - - /* local declarations */ - int i,j; - - /* node data: */ - const int numgrids=3; - const int NDOF1=1; - const int numdof=NDOF1*numgrids; - double xyz_list[numgrids][3]; - int doflist[numdof]; - int numberofdofspernode; - - /* gaussian points: */ - 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_l1l2l3[3]; - - /* matrices: */ - double L[numgrids]; - double B[2][numgrids]; - double Bprime[2][numgrids]; - double DL[2][2]={0.0}; - double DLprime[2][2]={0.0}; - double DL_scalar; - double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix - double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. - double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. - double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. - - double Jdettria; - - /*input parameters for structural analysis (diagnostic): */ - double dvx[2]; - double dvy[2]; - double vx,vy; - double dvxdx,dvydy; - double v_gauss[2]={0.0}; - - - double K[2][2]={0.0}; - double KDL[2][2]={0.0}; - int dofs[2]={0,1}; - int found=0; - - /*parameters: */ - bool artdiff; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /*retrieve some parameters: */ - this->parameters->FindParam(&artdiff,ArtDiffEnum); - - //Create Artificial diffusivity once for all if requested - if(artdiff){ - //Get the Jacobian determinant - gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); - - //Build K matrix (artificial diffusivity matrix) - inputs->GetParameterAverage(&v_gauss[0],VxAverageEnum); - inputs->GetParameterAverage(&v_gauss[1],VyAverageEnum); - - K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); - K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); - } - - /* 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; igGetParameterValue(&vx, &gauss_l1l2l3[0],VxAverageEnum); - inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyAverageEnum); - - inputs->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0],VxAverageEnum); - inputs->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0],VyAverageEnum); - - dvxdx=dvx[0]; - dvydy=dvy[1]; - - DL_scalar=gauss_weight*Jdettria; - - //Create DL and DLprime matrix - DL[0][0]=DL_scalar*dvxdx; - DL[1][1]=DL_scalar*dvydy; - - DLprime[0][0]=DL_scalar*vx; - DLprime[1][1]=DL_scalar*vy; - - //Do the triple product tL*D*L. - //Ke_gg_thickness=B'*DLprime*Bprime; - - TripleMultiply( &B[0][0],2,numdof,1, - &DL[0][0],2,2,0, - &B[0][0],2,numdof,0, - &Ke_gg_thickness1[0][0],0); - - TripleMultiply( &B[0][0],2,numdof,1, - &DLprime[0][0],2,2,0, - &Bprime[0][0],2,numdof,0, - &Ke_gg_thickness2[0][0],0); - - /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ - for( i=0; iGetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); - inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); - - DL_scalar=-gauss_weight*Jdettria; - - DLprime[0][0]=DL_scalar*vx; - DLprime[1][1]=DL_scalar*vy; - - //Do the triple product tL*D*L. - TripleMultiply( &B[0][0],2,numdof,1, - &DLprime[0][0],2,2,0, - &Bprime[0][0],2,numdof,0, - &Ke_gg2[0][0],0); - - /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ - for( i=0; iparameters->FindParam(&artdiff,ArtDiffEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /*Modify z so that it reflects the surface*/ - inputs->GetParameterValues(&surface_list[0],&gaussgrids[0][0],3,SurfaceEnum); - for(i=0;iGetParameterAverage(&nx,VxAverageEnum); - inputs->GetParameterAverage(&ny,VyAverageEnum); - if(nx==0 && ny==0){ - SurfaceNormal(&surface_normal[0],xyz_list); - nx=surface_normal[0]; - ny=surface_normal[1]; - } - if(nx==0 && ny==0){ - nx=0; - ny=1; - } - norm=pow( pow(nx,2)+pow(ny,2) , (double).5); - nx=nx/norm; - ny=ny/norm; - - //Create Artificial diffusivity once for all if requested - if(artdiff){ - //Get the Jacobian determinant - gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); - - //Build K matrix (artificial diffusivity matrix) - inputs->GetParameterAverage(&v_gauss[0],VxAverageEnum); - inputs->GetParameterAverage(&v_gauss[1],VyAverageEnum); - - K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!! - K[1][1]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); - } - - /* 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; igGetParameterValue(&vx,&gauss_l1l2l3[0],VxAverageEnum); - inputs->GetParameterValue(&vy,&gauss_l1l2l3[0],VyAverageEnum); - - inputs->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0],VxAverageEnum); - inputs->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0],VyAverageEnum); - - dvxdx=dvx[0]; - dvydy=dvy[1]; - - DL_scalar=gauss_weight*Jdettria; - - DLprime[0][0]=DL_scalar*nx; - DLprime[1][1]=DL_scalar*ny; - - //Do the triple product tL*D*L. - //Ke_gg_velocities=B'*DLprime*Bprime; - TripleMultiply( &B[0][0],2,numdof,1, - &DLprime[0][0],2,2,0, - &Bprime[0][0],2,numdof,0, - &Ke_gg_velocities2[0][0],0); - - /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ - for( i=0; iGetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); - - /*retrieve some parameters: */ - this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); - - /*First, if we are on water, return empty matrix: */ - if(onwater) return; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* 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; igGetParameterValue(&thickness, gauss_l1l2l3,ThicknessEnum); - - /*Get strain rate from velocity: */ - inputs->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss_l1l2l3,VxEnum,VyEnum); - inputs->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss_l1l2l3,VxOldEnum,VyOldEnum); - - /*Get viscosity: */ - matice->GetViscosity2d(&viscosity, &epsilon[0]); - matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant - onto this scalar matrix, so that we win some computational time: */ - newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - D_scalar=newviscosity*thickness*gauss_weight*Jdet; - - for (i=0;i<3;i++){ - D[i][i]=D_scalar; - } - - /*Get B and Bprime matrices: */ - GetB(&B[0][0], &xyz_list[0][0], gauss_l1l2l3); - GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); - - /* Do the triple product tB*D*Bprime: */ - TripleMultiply( &B[0][0],3,numdof,1, - &D[0][0],3,3,0, - &Bprime[0][0],3,numdof,0, - &Ke_gg_gaussian[0][0],0); - - /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ - for( i=0; iFindParam(&analysis_type,AnalysisTypeEnum); - - /*retrieve inputs :*/ - inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); - inputs->GetParameterValue(&drag_type,DragTypeEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - if (shelf){ - /*no friction, do nothing*/ - return; - } - - /*build friction object, used later on: */ - if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); - friction=new Friction("2d",inputs,matpar,analysis_type); - - /* 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; igGetAlpha2(&alpha2, gauss_l1l2l3,VxEnum,VyEnum,VzEnum); - - // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case, - //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ - inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],SurfaceEnum); - slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); - - if (slope_magnitude>MAXSLOPE){ - alpha2=pow((double)10,MOUNTAINKEXPONENT); - } - - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /*Get L matrix: */ - GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); - - - DL_scalar=alpha2*gauss_weight*Jdet; - for (i=0;i<2;i++){ - DL[i][i]=DL_scalar; - } - - /* Do the triple producte tL*D*L: */ - TripleMultiply( &L[0][0],2,numdof,1, - &DL[0][0],2,2,0, - &L[0][0],2,numdof,0, - &Ke_gg_gaussian[0][0],0); - - for( i=0; iGetParameterValue(&onwater,ElementOnWaterEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*Spawn 3 sing elements: */ - for(i=0;i<3;i++){ - sing=(Sing*)SpawnSing(i); - sing->CreateKMatrix(Kgg); - } - - /*clean up*/ - delete sing; - -} -/*}}}*/ -/*FUNCTION Tria::CreateKMatrixDiagnosticSurfaceVert {{{1*/ -void Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg){ - - int i,j; - - /* node data: */ - const int numgrids=3; - const int NDOF1=1; - const int numdof=NDOF1*numgrids; - double xyz_list[numgrids][3]; - int doflist[numdof]; - int numberofdofspernode; - - /* gaussian points: */ - 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_l1l2l3[3]; - - - /* surface normal: */ - double x4,y4,z4; - double x5,y5,z5; - double x6,y6,z6; - double v46[3]; - double v56[3]; - double normal[3]; - double norm_normal; - double nz; - - /*Matrices: */ - double DL_scalar; - double L[3]; - double Jdet; - - /* local element matrices: */ - double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix - double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* 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); - - /*Build normal vector to the surface:*/ - - x4=xyz_list[0][0]; - y4=xyz_list[0][1]; - z4=xyz_list[0][2]; - - x5=xyz_list[1][0]; - y5=xyz_list[1][1]; - z5=xyz_list[1][2]; - - x6=xyz_list[2][0]; - y6=xyz_list[2][1]; - z6=xyz_list[2][2]; - - v46[0]=x4-x6; - v46[1]=y4-y6; - v46[2]=z4-z6; - - v56[0]=x5-x6; - v56[1]=y5-y6; - v56[2]=z5-z6; - - normal[0]=(y4-y6)*(z5-z6)-(z4-z6)*(y5-y6); - normal[1]=(z4-z6)*(x5-x6)-(x4-x6)*(z5-z6); - normal[2]=(x4-x6)*(y5-y6)-(y4-y6)*(x5-x6); - - norm_normal=sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2)); - nz=1.0/norm_normal*normal[2]; - - /* Start looping on the number of gaussian points: */ - for (ig=0; igGetLatentHeat(); - heatcapacity=matpar->GetHeatCapacity(); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* Get gaussian points and weights: */ - GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - /* Start looping on the number of gauss (nodes on the bedrock) */ - for (ig=0; igparameters->FindParam(&dt,DtEnum); - this->parameters->FindParam(&artdiff,ArtDiffEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - //Create Artificial diffusivity once for all if requested - if(artdiff){ - //Get the Jacobian determinant - gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; - GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); - - //Build K matrix (artificial diffusivity matrix) - inputs->GetParameterAverage(&v_gauss[0],VxAverageEnum); - inputs->GetParameterAverage(&v_gauss[1],VyAverageEnum); - - K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); - K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); - } - - /* 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; igGetParameterValue(&vx,&gauss_l1l2l3[0],VxAverageEnum); - inputs->GetParameterValue(&vy,&gauss_l1l2l3[0],VyAverageEnum); - - inputs->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0],VxAverageEnum); - inputs->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0],VyAverageEnum); - - dvxdx=dvx[0]; - dvydy=dvy[1]; - - DL_scalar=dt*gauss_weight*Jdettria; - - //Create DL and DLprime matrix - DL[0][0]=DL_scalar*dvxdx; - DL[1][1]=DL_scalar*dvydy; - - DLprime[0][0]=DL_scalar*vx; - DLprime[1][1]=DL_scalar*vy; - - //Do the triple product tL*D*L. - //Ke_gg_thickness=B'*DL*B+B'*DLprime*Bprime; - - TripleMultiply( &B[0][0],2,numdof,1, - &DL[0][0],2,2,0, - &B[0][0],2,numdof,0, - &Ke_gg_thickness1[0][0],0); - - TripleMultiply( &B[0][0],2,numdof,1, - &DLprime[0][0],2,2,0, - &Bprime[0][0],2,numdof,0, - &Ke_gg_thickness2[0][0],0); - - /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ - for( i=0; iparameters->FindParam(&dt,DtEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* 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; igGetParameterValue(&vx,&gauss_l1l2l3[0],VxAverageEnum); - inputs->GetParameterValue(&vy,&gauss_l1l2l3[0],VyAverageEnum); - - DL_scalar=-dt*gauss_weight*Jdettria; - - DLprime[0][0]=DL_scalar*vx; - DLprime[1][1]=DL_scalar*vy; - - //Do the triple product tL*D*L. - TripleMultiply( &B[0][0],2,numdof,1, - &DLprime[0][0],2,2,0, - &Bprime[0][0],2,numdof,0, - &Ke_gg2[0][0],0); - - /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ - for( i=0; iparameters->FindParam(&dt,DtEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - //recover material parameters - mixed_layer_capacity=matpar->GetMixedLayerCapacity(); - thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - - - GaussTria (&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - /* Start looping on the number of gauss (nodes on the bedrock) */ - for (ig=0; igGetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); - inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); - - /* Add value into pe_g: */ - for( i=0; iGetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); - inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); - inputs->GetParameterValue(&dhdt_g, &gauss_l1l2l3[0],DhDtEnum); - - /* Add value into pe_g: */ - for( i=0; iGetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); - inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); - - /* Add value into pe_g: */ - for( i=0; iGetParameterValue(&meltingvalue, &gauss_l1l2l3[0],MeltingRateEnum); - - /*Get velocity at gaussian point: */ - inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); - inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); - - /*Get bed slope: */ - inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],BedEnum); - dbdx=slope[0]; - dbdy=slope[1]; - - /* Get Jacobian determinant: */ - GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - //Get L matrix if viscous basal drag present: - GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,NDOF1); - - - /*Build gaussian vector: */ - for(i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); - inputs->GetParameterValue(&drag_type,DragTypeEnum); - - /*First, if we are on water, return empty vector: */ - if(onwater)return; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - - /* Get gaussian points and weights: */ - GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/ - - /* Start looping on the number of gaussian points: */ - for (ig=0; igGetParameterValue(&thickness, &gauss_l1l2l3[0],ThicknessEnum); - inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],SurfaceEnum); - - /*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the - * element itself: */ - if(drag_type==1){ - inputs->GetParameterValue(&plastic_stress, &gauss_l1l2l3[0],DragCoefficientEnum); - } - - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /*Get nodal functions: */ - GetNodalFunctions(l1l2l3, gauss_l1l2l3); - - /*Compute driving stress: */ - driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness; - - /*Build pe_g_gaussian vector: */ - if(drag_type==1){ - for (i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); - - /*If on water, skip: */ - if(onwater)return; - - /*Spawn 3 sing elements: */ - for(i=0;i<3;i++){ - sing=(Sing*)SpawnSing(i); - sing->CreatePVector(pg); - } - - /*clean up*/ - delete sing; - -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorPrognostic {{{1*/ -void Tria::CreatePVectorPrognostic(Vec pg){ - - - /* local declarations */ - int i,j; - - /* node data: */ - const int numgrids=3; - const int NDOF1=1; - const int numdof=NDOF1*numgrids; - double xyz_list[numgrids][3]; - int doflist[numdof]; - int numberofdofspernode; - - /* gaussian points: */ - 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_l1l2l3[3]; - - /* matrix */ - double pe_g[numgrids]={0.0}; - double L[numgrids]; - double Jdettria; - - /*input parameters for structural analysis (diagnostic): */ - double accumulation_g; - double melting_g; - double thickness_g; - - /*parameters: */ - double dt; - - /*retrieve some parameters: */ - this->parameters->FindParam(&dt,DtEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* 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; igGetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); - inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); - inputs->GetParameterValue(&thickness_g, &gauss_l1l2l3[0],ThicknessEnum); - - /* Add value into pe_g: */ - for( i=0; iparameters->FindParam(&dt,DtEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - /* 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; igGetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); - inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); - inputs->GetParameterValue(&thickness_g, &gauss_l1l2l3[0],ThicknessEnum); - - /* Add value into pe_g: */ - for( i=0; iFindParam(&sub_analysis_type,AnalysisTypeEnum); - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - - /* Get gaussian points and weights: */ - GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/ - - - /* Start looping on the number of gaussian points: */ - for (ig=0; igGetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],SurfaceEnum); - } - if ( (sub_analysis_type==BedSlopeXAnalysisEnum) || (sub_analysis_type==BedSlopeYAnalysisEnum)){ - inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],BedEnum); - } - - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /*Get nodal functions: */ - GetNodalFunctions(l1l2l3, gauss_l1l2l3); - - /*Build pe_g_gaussian vector: */ - if ( (sub_analysis_type==SurfaceSlopeXAnalysisEnum) || (sub_analysis_type==BedSlopeXAnalysisEnum)){ - for(i=0;iGetMixedLayerCapacity(); - thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); - rho_water=matpar->GetRhoWater(); - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - beta=matpar->GetBeta(); - meltingpoint=matpar->GetMeltingPoint(); - - /*retrieve some solution parameters: */ - this->parameters->FindParam(&dt,DtEnum); - - /* Ice/ocean heat exchange flux on ice shelf base */ - - GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - /* Start looping on the number of gauss 2d (nodes on the bedrock) */ - for (ig=0; igGetParameterValue(&pressure, &gauss_coord[0],PressureEnum); - t_pmp=meltingpoint-beta*pressure; - - /*Calculate scalar parameter*/ - scalar_ocean=gauss_weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice); - if(dt){ - scalar_ocean=dt*scalar_ocean; - } - - for(i=0;i<3;i++){ - P_terms[i]+=scalar_ocean*l1l2l3[i]; - } - } - - /*Add pe_g to global vector pg: */ - VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES); - - cleanup_and_return: - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); - -} -/*}}}*/ -/*FUNCTION Tria::CreatePVectorThermalSheet {{{1*/ -void Tria::CreatePVectorThermalSheet( Vec pg){ - - int i,found; - - const int numgrids=3; - const int NDOF1=1; - const int numdof=numgrids*NDOF1; - int doflist[numdof]; - int numberofdofspernode; - double xyz_list[numgrids][3]; - - double rho_ice; - double heatcapacity; - - /*inputs: */ - double dt; - double pressure_list[3]; - double pressure; - int drag_type; - double basalfriction; - Friction* friction=NULL; - double alpha2,vx,vy; - double geothermalflux_value; - - /* gaussian points: */ - int num_area_gauss,ig; - double* gauss_weights = NULL; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double gauss_weight; - double gauss_coord[3]; - - /*matrices: */ - double Jdet; - double P_terms[numdof]={0.0}; - double l1l2l3[numgrids]; - double scalar; - - int analysis_type; - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList(&doflist[0],&numberofdofspernode); - - //recover material parameters - rho_ice=matpar->GetRhoIce(); - heatcapacity=matpar->GetHeatCapacity(); - - /*retrieve some parameters: */ - this->parameters->FindParam(&dt,DtEnum); - - /*Build frictoin element, needed later: */ - inputs->GetParameterValue(&drag_type,DragTypeEnum); - if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); - friction=new Friction("3d",inputs,matpar,analysis_type); - - /* Ice/ocean heat exchange flux on ice shelf base */ - GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); - - /* Start looping on the number of gauss 2d (nodes on the bedrock) */ - for (ig=0; igGetParameterValue(&geothermalflux_value, &gauss_coord[0],GeothermalFluxEnum); - - friction->GetAlpha2(&alpha2,&gauss_coord[0],VxEnum,VyEnum,VzEnum); - inputs->GetParameterValue(&vx, &gauss_coord[0],VxEnum); - inputs->GetParameterValue(&vy, &gauss_coord[0],VyEnum); - basalfriction= alpha2*(pow(vx,(double)2.0)+pow(vy,(double)2.0)); - - /*Calculate scalar parameter*/ - scalar=gauss_weight*Jdet*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice); - if(dt){ - scalar=dt*scalar; - } - - for(i=0;i<3;i++){ - P_terms[i]+=scalar*l1l2l3[i]; - } - } - - /*Add pe_g to global vector pg: */ - VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES); - - cleanup_and_return: - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); - delete friction; } @@ -3674,120 +909,4 @@ } /*}}}*/ -/*FUNCTION Tria::GetArea {{{1*/ -double Tria::GetArea(void){ - - double area=0; - const int numgrids=3; - double xyz_list[numgrids][3]; - double x1,y1,x2,y2,x3,y3; - - /*Get xyz list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - x1=xyz_list[0][0]; y1=xyz_list[0][1]; - x2=xyz_list[1][0]; y2=xyz_list[1][1]; - x3=xyz_list[2][0]; y3=xyz_list[2][1]; - - return x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1; -} -/*}}}*/ -/*FUNCTION Tria::GetAreaCoordinate {{{1*/ -double Tria::GetAreaCoordinate(double x, double y, int which_one){ - - double area=0; - const int numgrids=3; - double xyz_list[numgrids][3]; - double x1,y1,x2,y2,x3,y3; - - /*Get area: */ - area=this->GetArea(); - - /*Get xyz list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - x1=xyz_list[0][0]; y1=xyz_list[0][1]; - x2=xyz_list[1][0]; y2=xyz_list[1][1]; - x3=xyz_list[2][0]; y3=xyz_list[2][1]; - - if(which_one==1){ - /*Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area*/ - return ((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area; - } - else if(which_one==2){ - /*Get second area coordinate = det(x1-x3 x-x3 ; y1-y3 y-y3)/area*/ - return ((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area; - } - else if(which_one==3){ - /*Get third area coordinate 1-area1-area2: */ - return 1-((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area -((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area; - } - else ISSMERROR("%s%i%s\n"," error message: area coordinate ",which_one," done not exist!"); -} -/*}}}*/ -/*FUNCTION Tria::GetB {{{1*/ - -void Tria::GetB(double* B, double* xyz_list, double* gauss_l1l2l3){ - - /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. - * For grid i, Bi can be expressed in the actual coordinate system - * by: - * Bi=[ dh/dx 0 ] - * [ 0 dh/dy ] - * [ 1/2*dh/dy 1/2*dh/dx ] - * where h is the interpolation function for grid i. - * - * We assume B has been allocated already, of size: 3x(NDOF2*numgrids) - */ - - int i; - const int NDOF2=2; - const int numgrids=3; - - double dh1dh3[NDOF2][numgrids]; - - - /*Get dh1dh2dh3 in actual coordinate system: */ - GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss_l1l2l3); - - /*Build B: */ - for (i=0;iGetDofList(&doflist_per_node[0],&numberofdofspernode); - for(j=0;jGetDofList1(); - } - -} -/*}}}*/ -/*FUNCTION Tria::GetJacobian {{{1*/ -void Tria::GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3){ - - /*The Jacobian is constant over the element, discard the gaussian points. - * J is assumed to have been allocated of size NDOF2xNDOF2.*/ - - const int NDOF2=2; - const int numgrids=3; - double x1,y1,x2,y2,x3,y3; - - x1=*(xyz_list+numgrids*0+0); - y1=*(xyz_list+numgrids*0+1); - x2=*(xyz_list+numgrids*1+0); - y2=*(xyz_list+numgrids*1+1); - x3=*(xyz_list+numgrids*2+0); - y3=*(xyz_list+numgrids*2+1); - - - *(J+NDOF2*0+0)=0.5*(x2-x1); - *(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2); - *(J+NDOF2*0+1)=0.5*(y2-y1); - *(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2); -} -/*}}}*/ -/*FUNCTION Tria::GetJacobianDeterminant2d {{{1*/ -void Tria::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss_l1l2l3){ - - /*The Jacobian determinant is constant over the element, discard the gaussian points. - * J is assumed to have been allocated of size NDOF2xNDOF2.*/ - - double x1,x2,x3,y1,y2,y3; - - x1=*(xyz_list+3*0+0); - y1=*(xyz_list+3*0+1); - x2=*(xyz_list+3*1+0); - y2=*(xyz_list+3*1+1); - x3=*(xyz_list+3*2+0); - y3=*(xyz_list+3*2+1); - - - *Jdet=SQRT3/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)); - - - if(Jdet<0){ - ISSMERROR("negative jacobian determinant!"); - } - -} -/*}}}*/ -/*FUNCTION Tria::GetJacobianDeterminant3d {{{1*/ -void Tria::GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss_l1l2l3){ - - /*The Jacobian determinant is constant over the element, discard the gaussian points. - * J is assumed to have been allocated of size NDOF2xNDOF2.*/ - - double x1,x2,x3,y1,y2,y3,z1,z2,z3; - - x1=*(xyz_list+3*0+0); - y1=*(xyz_list+3*0+1); - z1=*(xyz_list+3*0+2); - x2=*(xyz_list+3*1+0); - y2=*(xyz_list+3*1+1); - z2=*(xyz_list+3*1+2); - x3=*(xyz_list+3*2+0); - y3=*(xyz_list+3*2+1); - z3=*(xyz_list+3*2+2); - - - *Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5); - - if(Jdet<0){ - ISSMERROR("negative jacobian determinant!"); - } - -} -/*}}}*/ -/*FUNCTION Tria::GetJacobianInvert {{{1*/ -void Tria::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss_l1l2l3){ - - double Jdet; - const int NDOF2=2; - const int numgrids=3; - - /*Call Jacobian routine to get the jacobian:*/ - GetJacobian(Jinv, xyz_list, gauss_l1l2l3); - - /*Invert Jacobian matrix: */ - MatrixInverse(Jinv,NDOF2,NDOF2,NULL,0,&Jdet); - -} -/*}}}*/ -/*FUNCTION Tria::GetL {{{1*/ - -void Tria::GetL(double* L, double* xyz_list, double* gauss_l1l2l3,int numdof){ - - /*Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof. - * For grid i, Li can be expressed in the actual coordinate system - * by: - * numdof=1: - * Li=h; - * numdof=2: - * Li=[ h 0 ] - * [ 0 h ] - * where h is the interpolation function for grid i. - * - * We assume L has been allocated already, of size: numgrids (numdof=1), or numdofx(numdof*numgrids) (numdof=2) - */ - - int i; - const int NDOF2=2; - const int numgrids=3; - - double l1l2l3[3]; - - - /*Get l1l2l3 in actual coordinate system: */ - GetNodalFunctions(l1l2l3, gauss_l1l2l3); - -#ifdef _DELUG_ - for (i=0;i<3;i++){ - printf("Node %i h=%lf \n",i,l1l2l3[i]); - } -#endif - - /*Build L: */ - if(numdof==1){ - 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==DiagnosticHorizAnalysisEnum) + GetSolutionFromInputsDiagnosticHoriz(solution); + else + ISSMERROR("%s%i%s\n","analysis: ",analysis_type," not supported yet"); + +} +/*}}}*/ /*FUNCTION Tria::GetThicknessList {{{1*/ void Tria::GetThicknessList(double* thickness_list){ @@ -4207,4 +983,23 @@ } /*}}}*/ +/*FUNCTION Tria::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/ +void Tria::GetVectorFromInputs(Vec vector,int NameEnum){ + + int i; + const int numvertices=3; + 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 Tria::Gradj {{{1*/ void Tria::Gradj(Vec gradient,int control_type){ @@ -4226,339 +1021,4 @@ } else ISSMERROR("%s%i","control type not supported yet: ",control_type); -} -/*}}}*/ -/*FUNCTION Tria::GradjDrag {{{1*/ -void Tria::GradjDrag(Vec gradient){ - - - int i; - - /* node data: */ - const int numgrids=3; - const int NDOF2=2; - const int numdof=NDOF2*numgrids; - double xyz_list[numgrids][3]; - int doflist1[numgrids]; - double dh1dh3[NDOF2][numgrids]; - - /* grid data: */ - double adjx_list[numgrids]; - double adjy_list[numgrids]; - - /* gaussian points: */ - 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_l1l2l3[3]; - - /* parameters: */ - double dk[NDOF2]; - double vx,vy; - double lambda,mu; - double bed,thickness,Neff; - double alpha_complement; - int drag_type; - double drag; - Friction* friction=NULL; - - /*element vector at the gaussian points: */ - double grade_g[numgrids]={0.0}; - double grade_g_gaussian[numgrids]; - - /* Jacobian: */ - double Jdet; - - /*nodal functions: */ - double l1l2l3[3]; - - /* strain rate: */ - double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - - /*inputs: */ - bool shelf; - - /*parameters: */ - double cm_noisedmp; - double cm_mindmp_slope; - double cm_mindmp_value; - double cm_maxdmp_value; - double cm_maxdmp_slope; - - int analysis_type; - - /*retrive parameters: */ - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*retrieve inputs :*/ - inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); - - /*retrieve some parameters: */ - this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); - this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum); - this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum); - this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum); - this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum); - - - /*Get out if shelf*/ - if(shelf)return; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList1(&doflist1[0]); - - /*Build frictoin element, needed later: */ - inputs->GetParameterValue(&drag_type,DragTypeEnum); - friction=new Friction("2d",inputs,matpar,analysis_type); - - /* 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, 4); - - /* Start looping on the number of gaussian points: */ - for (ig=0; igGetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum); - else alpha_complement=0; - - /*Recover alpha_complement and k: */ - inputs->GetParameterValue(&drag, gauss_l1l2l3,DragCoefficientEnum); - - /*recover lambda and mu: */ - inputs->GetParameterValue(&lambda, gauss_l1l2l3,AdjointxEnum); - inputs->GetParameterValue(&mu, gauss_l1l2l3,AdjointyEnum); - - /*recover vx and vy: */ - inputs->GetParameterValue(&vx, gauss_l1l2l3,VxEnum); - inputs->GetParameterValue(&vy, gauss_l1l2l3,VyEnum); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /* Get nodal functions value at gaussian point:*/ - GetNodalFunctions(l1l2l3, gauss_l1l2l3); - - /*Get nodal functions derivatives*/ - GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3); - - /*Get k derivative: dk/dx */ - inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); - - /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for (i=0;icm_maxdmp_value){ - grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; - } - } - - /*Add gradje_g_gaussian vector to gradje_g: */ - for( i=0; iFindParam(&analysis_type,AnalysisTypeEnum); - - /*retrieve inputs :*/ - inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); - inputs->GetParameterValue(&drag_type,DragTypeEnum); - - /*retrieve some parameters: */ - this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); - this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum); - this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum); - this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum); - this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum); - - /*Get out if shelf*/ - if(shelf)return; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - GetDofList1(&doflist1[0]); - - /*Build frictoin element, needed later: */ - inputs->GetParameterValue(&drag_type,DragTypeEnum); - friction=new Friction("2d",inputs,matpar,analysis_type); - - /* 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, 4); - - /* Start looping on the number of gaussian points: */ - for (ig=0; igGetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum); - else alpha_complement=0; - inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum); - - /*recover lambda mu and xi: */ - inputs->GetParameterValue(&lambda, &gauss_l1l2l3[0],AdjointxEnum); - inputs->GetParameterValue(&mu, &gauss_l1l2l3[0],AdjointyEnum); - inputs->GetParameterValue(&xi, &gauss_l1l2l3[0],AdjointzEnum); - - /*recover vx vy and vz: */ - inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); - inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); - inputs->GetParameterValue(&vz, &gauss_l1l2l3[0],VzEnum); - - /*Get normal vector to the bed */ - SurfaceNormal(&surface_normal[0],xyz_list); - - 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]; - - /* Get Jacobian determinant: */ - GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); - - /* Get nodal functions value at gaussian point:*/ - GetNodalFunctions(l1l2l3, gauss_l1l2l3); - - /*Get nodal functions derivatives*/ - GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3); - - /*Get k derivative: dk/dx */ - inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); - - /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ - for (i=0;icm_maxdmp_value){ - grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; - } - } - - /*Add gradje_g_gaussian vector to gradje_g: */ - for( i=0; iinputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0])); - - cleanup_and_return: - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&gauss_weights); - delete friction; - } /*}}}*/ @@ -4709,4 +1169,299 @@ } /*}}}*/ +/*FUNCTION Tria::GradjDrag {{{1*/ +void Tria::GradjDrag(Vec gradient){ + + + int i; + + /* node data: */ + const int numgrids=3; + const int NDOF2=2; + const int numdof=NDOF2*numgrids; + double xyz_list[numgrids][3]; + int doflist1[numgrids]; + double dh1dh3[NDOF2][numgrids]; + + /* grid data: */ + double adjx_list[numgrids]; + double adjy_list[numgrids]; + + /* gaussian points: */ + 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_l1l2l3[3]; + + /* parameters: */ + double dk[NDOF2]; + double vx,vy; + double lambda,mu; + double bed,thickness,Neff; + double alpha_complement; + int drag_type; + double drag; + Friction* friction=NULL; + + /*element vector at the gaussian points: */ + double grade_g[numgrids]={0.0}; + double grade_g_gaussian[numgrids]; + + /* Jacobian: */ + double Jdet; + + /*nodal functions: */ + double l1l2l3[3]; + + /* strain rate: */ + double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + + /*inputs: */ + bool shelf; + + /*parameters: */ + double cm_noisedmp; + double cm_mindmp_slope; + double cm_mindmp_value; + double cm_maxdmp_value; + double cm_maxdmp_slope; + + int analysis_type; + + /*retrive parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*retrieve inputs :*/ + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); + this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum); + this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum); + this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum); + this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum); + + + /*Get out if shelf*/ + if(shelf)return; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList1(&doflist1[0]); + + /*Build frictoin element, needed later: */ + inputs->GetParameterValue(&drag_type,DragTypeEnum); + friction=new Friction("2d",inputs,matpar,analysis_type); + + /* 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, 4); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igGetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum); + else alpha_complement=0; + + /*Recover alpha_complement and k: */ + inputs->GetParameterValue(&drag, gauss_l1l2l3,DragCoefficientEnum); + + /*recover lambda and mu: */ + inputs->GetParameterValue(&lambda, gauss_l1l2l3,AdjointxEnum); + inputs->GetParameterValue(&mu, gauss_l1l2l3,AdjointyEnum); + + /*recover vx and vy: */ + inputs->GetParameterValue(&vx, gauss_l1l2l3,VxEnum); + inputs->GetParameterValue(&vy, gauss_l1l2l3,VyEnum); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /* Get nodal functions value at gaussian point:*/ + GetNodalFunctions(l1l2l3, gauss_l1l2l3); + + /*Get nodal functions derivatives*/ + GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3); + + /*Get k derivative: dk/dx */ + inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for (i=0;icm_maxdmp_value){ + grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; + } + } + + /*Add gradje_g_gaussian vector to gradje_g: */ + for( i=0; iinputs->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!"); + + /*Scale: */ + yinput->AXPY(xinput,scalar); +} +/*}}}*/ +/*FUNCTION Tria::InputControlConstrain(int control_type, double cm_min, double cm_max){{{1*/ +void Tria::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 Tria::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/ +void Tria::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 Tria::InputDepthAverageAtBase {{{1*/ +void Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){ + + /*New input*/ + Input* oldinput=NULL; + Input* newinput=NULL; + + /*copy input of enum_type*/ + oldinput=this->inputs->GetInput(enum_type); + if(!oldinput)ISSMERROR("%s%s"," could not find old input with enum: ",EnumAsString(enum_type)); + newinput=(Input*)oldinput->copy(); + + /*Assign new name (average)*/ + newinput->ChangeEnum(average_enum_type); + + /*Add new input to current element*/ + this->inputs->AddInput(newinput); + +} +/*}}}*/ +/*FUNCTION Tria::InputDuplicate(int original_enum,int new_enum){{{1*/ +void Tria::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 Tria::InputScale(int enum_type,double scale_factor){{{1*/ +void Tria::InputScale(int enum_type,double scale_factor){ + + Input* input=NULL; + + /*Make a copy of the original input: */ + input=(Input*)this->inputs->GetInput(enum_type); + + /*Scale: */ + input->Scale(scale_factor); +} +/*}}}*/ +/*FUNCTION Tria::InputToResult(int enum_type,int step,double time){{{1*/ +void Tria::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 Tria::MassFlux {{{1*/ double Tria::MassFlux( double* segment){ @@ -4766,4 +1521,324 @@ ); return mass_flux; +} +/*}}}*/ +/*FUNCTION Tria::MaxAbsVx(double* pmaxabsvx, bool process_units);{{{1*/ +void Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){ + + int i; + int dim; + const int numgrids=3; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,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 Tria::MaxAbsVy(double* pmaxabsvy, bool process_units);{{{1*/ +void Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){ + + int i; + int dim; + const int numgrids=3; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,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 Tria::MaxAbsVz(double* pmaxabsvz, bool process_units);{{{1*/ +void Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){ + + int i; + int dim; + const int numgrids=3; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,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 Tria::MaxVel(double* pmaxvel, bool process_units);{{{1*/ +void Tria::MaxVel(double* pmaxvel, bool process_units){ + + int i; + int dim; + const int numgrids=3; + double gaussgrids[numgrids][3]={{1,0,0},{0,1,0},{0,0,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 Tria::MaxVx(double* pmaxvx, bool process_units);{{{1*/ +void Tria::MaxVx(double* pmaxvx, bool process_units){ + + int i; + int dim; + const int numgrids=3; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,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 Tria::MaxVy(double* pmaxvy, bool process_units);{{{1*/ +void Tria::MaxVy(double* pmaxvy, bool process_units){ + + int i; + int dim; + const int numgrids=3; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,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 Tria::MaxVz(double* pmaxvz, bool process_units);{{{1*/ +void Tria::MaxVz(double* pmaxvz, bool process_units){ + + int i; + int dim; + const int numgrids=3; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,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 Tria::MinVel(double* pminvel, bool process_units);{{{1*/ +void Tria::MinVel(double* pminvel, bool process_units){ + + int i; + int dim; + const int numgrids=3; + double gaussgrids[numgrids][3]={{1,0,0},{0,1,0},{0,0,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(&fit,FitEnum); - inputs->GetParameterValue(&onwater,ElementOnWaterEnum); - - /*If fit!=3, do not compute surface: */ - if(fit!=3)return 0; - - /*If on water, return 0: */ - if(onwater)return 0; - - /* Get node coordinates and dof list: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); - - for (i=0;i<3;i++){ - v13[i]=xyz_list[0][i]-xyz_list[2][i]; - v23[i]=xyz_list[1][i]-xyz_list[2][i]; - } - - normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; - normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; - normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; - - S = 0.5 * sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2)); - - /*Return: */ - return S; -} -/*}}}*/ -/*FUNCTION Tria::InputUpdateFromVector(double* vector, int name, int type);{{{1*/ -void Tria::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[3]; - - /*Get values on the 6 vertices*/ - for (int i=0;i<3;i++){ - values[i]=vector[this->nodes[i]->GetVertexDof()]; - } - - /*update input*/ - this->inputs->AddInput(new TriaVertexInput(name,values)); - return; - - default: - - ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type)); - } -} -/*}}}*/ -/*FUNCTION Tria::InputUpdateFromVector(int* vector, int name, int type);{{{1*/ -void Tria::InputUpdateFromVector(int* vector, int name, int type){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Tria::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/ -void Tria::InputUpdateFromVector(bool* vector, int name, int type){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ -/*FUNCTION Tria::InputUpdateFromConstant(int value, int name);{{{1*/ -void Tria::InputUpdateFromConstant(int constant, int name){ - /*Nothing updated for now*/ -} -/*}}}*/ -/*FUNCTION Tria::InputUpdateFromConstant(double value, int name);{{{1*/ -void Tria::InputUpdateFromConstant(double constant, int name){ - /*Nothing updated for now*/ -} -/*}}}*/ -/*FUNCTION Tria::InputUpdateFromConstant(bool value, int name);{{{1*/ -void Tria::InputUpdateFromConstant(bool constant, int name){ - /*Nothing updated for now*/ -} -/*}}}*/ -/*FUNCTION Tria::InputDepthAverageAtBase {{{1*/ -void Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){ - - /*New input*/ - Input* oldinput=NULL; - Input* newinput=NULL; - - /*copy input of enum_type*/ - oldinput=this->inputs->GetInput(enum_type); - if(!oldinput)ISSMERROR("%s%s"," could not find old input with enum: ",EnumAsString(enum_type)); - newinput=(Input*)oldinput->copy(); - - /*Assign new name (average)*/ - newinput->ChangeEnum(average_enum_type); - - /*Add new input to current element*/ - this->inputs->AddInput(newinput); - + int row; + int vertices_ids[3]; + + + /*recover pointer: */ + row=*prow; + + /*will be needed later: */ + for(i=0;i<3;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch. + + for(i=0;iresults->Size();i++){ + ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); + + /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand + *it to the result object, to fill the rest: */ + patch->fillelementinfo(row,this->id,vertices_ids,3); + elementresult->PatchFill(row,patch); + + /*increment rower: */ + row++; + } + + /*Assign output pointers:*/ + *prow=row; } /*}}}*/ @@ -5140,461 +2092,3313 @@ } /*}}}*/ -/*FUNCTION Tria::PatchFill(int* prow, Patch* patch){{{1*/ -void Tria::PatchFill(int* prow, Patch* patch){ +/*FUNCTION Tria::ProcessResultsUnits(void){{{1*/ +void Tria::ProcessResultsUnits(void){ int i; - int row; - int vertices_ids[3]; - - - /*recover pointer: */ - row=*prow; - - /*will be needed later: */ - for(i=0;i<3;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch. for(i=0;iresults->Size();i++){ ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i); - - /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand - *it to the result object, to fill the rest: */ - patch->fillelementinfo(row,this->id,vertices_ids,3); - elementresult->PatchFill(row,patch); - - /*increment rower: */ - row++; + elementresult->ProcessUnits(this->parameters); + } + +} +/*}}}*/ +/*FUNCTION Tria::SurfaceArea {{{1*/ +double Tria::SurfaceArea(void){ + + int i; + + /* output: */ + double S; + + /* node data: */ + int numgrids=3; + double xyz_list[numgrids][3]; + double v13[3]; + double v23[3]; + double normal[3]; + + /*inputs: */ + bool onwater; + int fit; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&fit,FitEnum); + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + + /*If fit!=3, do not compute surface: */ + if(fit!=3)return 0; + + /*If on water, return 0: */ + if(onwater)return 0; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + + for (i=0;i<3;i++){ + v13[i]=xyz_list[0][i]-xyz_list[2][i]; + v23[i]=xyz_list[1][i]-xyz_list[2][i]; + } + + normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; + normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; + normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; + + S = 0.5 * sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2)); + + /*Return: */ + return S; +} +/*}}}*/ +/*FUNCTION Tria::Update(IoModel* iomodel,int analysis_counter,int analysis_type){{{1*/ +void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index + + /*Intermediaries*/ + int i; + int tria_node_ids[3]; + int tria_vertex_ids[3]; + double nodeinputs[3]; + + /*Checks if debuging*/ + /*{{{2*/ + ISSMASSERT(iomodel->elements); + /*}}}*/ + + /*Recover vertices ids needed to initialize inputs*/ + for(i=0;i<3;i++){ + tria_vertex_ids[i]=(int)iomodel->elements[3*index+i]; //ids for vertices are in the elements array from Matlab + } + + /*Recover nodes ids needed to initialize the node hook.*/ + if (analysis_type==Prognostic2AnalysisEnum || analysis_type==Balancedthickness2AnalysisEnum){ + /*Discontinuous Galerkin*/ + tria_node_ids[0]=iomodel->nodecounter+3*index+1; + tria_node_ids[1]=iomodel->nodecounter+3*index+2; + tria_node_ids[2]=iomodel->nodecounter+3*index+3; + } + else{ + /*Continuous Galerkin*/ + for(i=0;i<3;i++){ + tria_node_ids[i]=iomodel->nodecounter+(int)*(iomodel->elements+3*index+i); //ids for vertices are in the elements array from Matlab + } + } + + /*hooks: */ + this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type + + /*add as many inputs per element as requested:*/ + if (iomodel->thickness) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,nodeinputs)); + } + if (iomodel->surface) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,nodeinputs)); + } + if (iomodel->bed) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(BedEnum,nodeinputs)); + } + if (iomodel->drag_coefficient) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->drag_coefficient[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(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<3;i++)nodeinputs[i]=iomodel->melting_rate[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(MeltingRateEnum,nodeinputs)); + } + if (iomodel->accumulation_rate) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(AccumulationRateEnum,nodeinputs)); + } + if (iomodel->geothermalflux) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->geothermalflux[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(GeothermalFluxEnum,nodeinputs)); + } + if (iomodel->dhdt) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->dhdt[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(DhDtEnum,nodeinputs)); + } + if (iomodel->pressure) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->pressure[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(PressureEnum,nodeinputs)); + } + if (iomodel->temperature) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->temperature[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(TemperatureEnum,nodeinputs)); + } + /*vx,vy and vz: */ + if (iomodel->vx) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(VxOldEnum,nodeinputs)); + } + if (iomodel->vy) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(VyOldEnum,nodeinputs)); + } + if (iomodel->vz) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(VzOldEnum,nodeinputs)); + } + if (iomodel->vx_obs) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(VxObsEnum,nodeinputs)); + } + if (iomodel->vy_obs) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(VyObsEnum,nodeinputs)); + } + if (iomodel->vz_obs) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(VzObsEnum,nodeinputs)); + } + if (iomodel->weights) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->weights[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(WeightsEnum,nodeinputs)); + } + if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index])); + if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index])); + if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index])); + if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index])); + + /*Defaults if not provided in iomodel*/ + switch(analysis_type){ + + case DiagnosticHorizAnalysisEnum: case DiagnosticVertAnalysisEnum: case DiagnosticStokesAnalysisEnum: + + /*default vx,vy and vz: either observation or 0 */ + if(!iomodel->vx){ + if (iomodel->vx_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts; + else for(i=0;i<3;i++)nodeinputs[i]=0; + this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(VxOldEnum,nodeinputs)); + } + if(!iomodel->vy){ + if (iomodel->vy_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts; + else for(i=0;i<3;i++)nodeinputs[i]=0; + this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(VyOldEnum,nodeinputs)); + } + if(!iomodel->vz){ + if (iomodel->vz_obs) for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_vertex_ids[i]-1]/iomodel->yts; + else for(i=0;i<3;i++)nodeinputs[i]=0; + this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(VzOldEnum,nodeinputs)); + } + break; + + default: + /*No update for other solution types*/ + break; + + } + + //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. + this->parameters=NULL; + +} +/*}}}*/ +/*FUNCTION Tria::UpdateGeometry{{{1*/ +void Tria::UpdateGeometry(void){ + + /*Intermediaries*/ + 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 + } + +} +/*}}}*/ + + +/*Tria specific routines: */ +/*FUNCTION Tria::CreateKMatrixBalancedthickness {{{1*/ +void Tria::CreateKMatrixBalancedthickness(Mat Kgg){ + + /* local declarations */ + int i,j; + + /* node data: */ + const int numgrids=3; + const int NDOF1=1; + const int numdof=NDOF1*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdof]; + int numberofdofspernode; + + /* gaussian points: */ + 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_l1l2l3[3]; + + /* matrices: */ + double L[numgrids]; + double B[2][numgrids]; + double Bprime[2][numgrids]; + double DL[2][2]={0.0}; + double DLprime[2][2]={0.0}; + double DL_scalar; + double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix + double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. + double Ke_gg_thickness1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. + double Ke_gg_thickness2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point. + + double Jdettria; + + /*input parameters for structural analysis (diagnostic): */ + double dvx[2]; + double dvy[2]; + double vx,vy; + double dvxdx,dvydy; + double v_gauss[2]={0.0}; + + + double K[2][2]={0.0}; + double KDL[2][2]={0.0}; + int dofs[2]={0,1}; + int found=0; + + /*parameters: */ + bool artdiff; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /*retrieve some parameters: */ + this->parameters->FindParam(&artdiff,ArtDiffEnum); + + //Create Artificial diffusivity once for all if requested + if(artdiff){ + //Get the Jacobian determinant + gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + + //Build K matrix (artificial diffusivity matrix) + inputs->GetParameterAverage(&v_gauss[0],VxAverageEnum); + inputs->GetParameterAverage(&v_gauss[1],VyAverageEnum); + + K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); + K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); + } + + /* 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; igGetParameterValue(&vx, &gauss_l1l2l3[0],VxAverageEnum); + inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyAverageEnum); + + inputs->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0],VxAverageEnum); + inputs->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0],VyAverageEnum); + + dvxdx=dvx[0]; + dvydy=dvy[1]; + + DL_scalar=gauss_weight*Jdettria; + + //Create DL and DLprime matrix + DL[0][0]=DL_scalar*dvxdx; + DL[1][1]=DL_scalar*dvydy; + + DLprime[0][0]=DL_scalar*vx; + DLprime[1][1]=DL_scalar*vy; + + //Do the triple product tL*D*L. + //Ke_gg_thickness=B'*DLprime*Bprime; + + TripleMultiply( &B[0][0],2,numdof,1, + &DL[0][0],2,2,0, + &B[0][0],2,numdof,0, + &Ke_gg_thickness1[0][0],0); + + TripleMultiply( &B[0][0],2,numdof,1, + &DLprime[0][0],2,2,0, + &Bprime[0][0],2,numdof,0, + &Ke_gg_thickness2[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ + for( i=0; iGetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); + inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); + + DL_scalar=-gauss_weight*Jdettria; + + DLprime[0][0]=DL_scalar*vx; + DLprime[1][1]=DL_scalar*vy; + + //Do the triple product tL*D*L. + TripleMultiply( &B[0][0],2,numdof,1, + &DLprime[0][0],2,2,0, + &Bprime[0][0],2,numdof,0, + &Ke_gg2[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ + for( i=0; iparameters->FindParam(&artdiff,ArtDiffEnum); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /*Modify z so that it reflects the surface*/ + inputs->GetParameterValues(&surface_list[0],&gaussgrids[0][0],3,SurfaceEnum); + for(i=0;iGetParameterAverage(&nx,VxAverageEnum); + inputs->GetParameterAverage(&ny,VyAverageEnum); + if(nx==0 && ny==0){ + SurfaceNormal(&surface_normal[0],xyz_list); + nx=surface_normal[0]; + ny=surface_normal[1]; + } + if(nx==0 && ny==0){ + nx=0; + ny=1; + } + norm=pow( pow(nx,2)+pow(ny,2) , (double).5); + nx=nx/norm; + ny=ny/norm; + + //Create Artificial diffusivity once for all if requested + if(artdiff){ + //Get the Jacobian determinant + gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + + //Build K matrix (artificial diffusivity matrix) + inputs->GetParameterAverage(&v_gauss[0],VxAverageEnum); + inputs->GetParameterAverage(&v_gauss[1],VyAverageEnum); + + K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!! + K[1][1]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); + } + + /* 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; igGetParameterValue(&vx,&gauss_l1l2l3[0],VxAverageEnum); + inputs->GetParameterValue(&vy,&gauss_l1l2l3[0],VyAverageEnum); + + inputs->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0],VxAverageEnum); + inputs->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0],VyAverageEnum); + + dvxdx=dvx[0]; + dvydy=dvy[1]; + + DL_scalar=gauss_weight*Jdettria; + + DLprime[0][0]=DL_scalar*nx; + DLprime[1][1]=DL_scalar*ny; + + //Do the triple product tL*D*L. + //Ke_gg_velocities=B'*DLprime*Bprime; + TripleMultiply( &B[0][0],2,numdof,1, + &DLprime[0][0],2,2,0, + &Bprime[0][0],2,numdof,0, + &Ke_gg_velocities2[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ + for( i=0; iGetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum); + + /*First, if we are on water, return empty matrix: */ + if(onwater) return; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* 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; igGetParameterValue(&thickness, gauss_l1l2l3,ThicknessEnum); + + /*Get strain rate from velocity: */ + inputs->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss_l1l2l3,VxEnum,VyEnum); + inputs->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss_l1l2l3,VxOldEnum,VyOldEnum); + + /*Get viscosity: */ + matice->GetViscosity2d(&viscosity, &epsilon[0]); + matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant + onto this scalar matrix, so that we win some computational time: */ + newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); + D_scalar=newviscosity*thickness*gauss_weight*Jdet; + + for (i=0;i<3;i++){ + D[i][i]=D_scalar; + } + + /*Get B and Bprime matrices: */ + GetB(&B[0][0], &xyz_list[0][0], gauss_l1l2l3); + GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); + + /* Do the triple product tB*D*Bprime: */ + TripleMultiply( &B[0][0],3,numdof,1, + &D[0][0],3,3,0, + &Bprime[0][0],3,numdof,0, + &Ke_gg_gaussian[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ + for( i=0; iFindParam(&analysis_type,AnalysisTypeEnum); + + /*retrieve inputs :*/ + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + inputs->GetParameterValue(&drag_type,DragTypeEnum); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + if (shelf){ + /*no friction, do nothing*/ + return; + } + + /*build friction object, used later on: */ + if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); + friction=new Friction("2d",inputs,matpar,analysis_type); + + /* 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; igGetAlpha2(&alpha2, gauss_l1l2l3,VxEnum,VyEnum,VzEnum); + + // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case, + //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ + inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],SurfaceEnum); + slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); + + if (slope_magnitude>MAXSLOPE){ + alpha2=pow((double)10,MOUNTAINKEXPONENT); + } + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /*Get L matrix: */ + GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); + + + DL_scalar=alpha2*gauss_weight*Jdet; + for (i=0;i<2;i++){ + DL[i][i]=DL_scalar; + } + + /* Do the triple producte tL*D*L: */ + TripleMultiply( &L[0][0],2,numdof,1, + &DL[0][0],2,2,0, + &L[0][0],2,numdof,0, + &Ke_gg_gaussian[0][0],0); + + for( i=0; iGetParameterValue(&onwater,ElementOnWaterEnum); + + /*If on water, skip: */ + if(onwater)return; + + /*Spawn 3 sing elements: */ + for(i=0;i<3;i++){ + sing=(Sing*)SpawnSing(i); + sing->CreateKMatrix(Kgg); + } + + /*clean up*/ + delete sing; + +} +/*}}}*/ +/*FUNCTION Tria::CreateKMatrixDiagnosticSurfaceVert {{{1*/ +void Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg){ + + int i,j; + + /* node data: */ + const int numgrids=3; + const int NDOF1=1; + const int numdof=NDOF1*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdof]; + int numberofdofspernode; + + /* gaussian points: */ + 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_l1l2l3[3]; + + + /* surface normal: */ + double x4,y4,z4; + double x5,y5,z5; + double x6,y6,z6; + double v46[3]; + double v56[3]; + double normal[3]; + double norm_normal; + double nz; + + /*Matrices: */ + double DL_scalar; + double L[3]; + double Jdet; + + /* local element matrices: */ + double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix + double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* 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); + + /*Build normal vector to the surface:*/ + + x4=xyz_list[0][0]; + y4=xyz_list[0][1]; + z4=xyz_list[0][2]; + + x5=xyz_list[1][0]; + y5=xyz_list[1][1]; + z5=xyz_list[1][2]; + + x6=xyz_list[2][0]; + y6=xyz_list[2][1]; + z6=xyz_list[2][2]; + + v46[0]=x4-x6; + v46[1]=y4-y6; + v46[2]=z4-z6; + + v56[0]=x5-x6; + v56[1]=y5-y6; + v56[2]=z5-z6; + + normal[0]=(y4-y6)*(z5-z6)-(z4-z6)*(y5-y6); + normal[1]=(z4-z6)*(x5-x6)-(x4-x6)*(z5-z6); + normal[2]=(x4-x6)*(y5-y6)-(y4-y6)*(x5-x6); + + norm_normal=sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2)); + nz=1.0/norm_normal*normal[2]; + + /* Start looping on the number of gaussian points: */ + for (ig=0; igGetLatentHeat(); + heatcapacity=matpar->GetHeatCapacity(); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Get gaussian points and weights: */ + GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gauss (nodes on the bedrock) */ + for (ig=0; igparameters->FindParam(&dt,DtEnum); + this->parameters->FindParam(&artdiff,ArtDiffEnum); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + //Create Artificial diffusivity once for all if requested + if(artdiff){ + //Get the Jacobian determinant + gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + + //Build K matrix (artificial diffusivity matrix) + inputs->GetParameterAverage(&v_gauss[0],VxAverageEnum); + inputs->GetParameterAverage(&v_gauss[1],VyAverageEnum); + + K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); + K[1][1]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); + } + + /* 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; igGetParameterValue(&vx,&gauss_l1l2l3[0],VxAverageEnum); + inputs->GetParameterValue(&vy,&gauss_l1l2l3[0],VyAverageEnum); + + inputs->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0],VxAverageEnum); + inputs->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0],VyAverageEnum); + + dvxdx=dvx[0]; + dvydy=dvy[1]; + + DL_scalar=dt*gauss_weight*Jdettria; + + //Create DL and DLprime matrix + DL[0][0]=DL_scalar*dvxdx; + DL[1][1]=DL_scalar*dvydy; + + DLprime[0][0]=DL_scalar*vx; + DLprime[1][1]=DL_scalar*vy; + + //Do the triple product tL*D*L. + //Ke_gg_thickness=B'*DL*B+B'*DLprime*Bprime; + + TripleMultiply( &B[0][0],2,numdof,1, + &DL[0][0],2,2,0, + &B[0][0],2,numdof,0, + &Ke_gg_thickness1[0][0],0); + + TripleMultiply( &B[0][0],2,numdof,1, + &DLprime[0][0],2,2,0, + &Bprime[0][0],2,numdof,0, + &Ke_gg_thickness2[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ + for( i=0; iparameters->FindParam(&dt,DtEnum); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* 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; igGetParameterValue(&vx,&gauss_l1l2l3[0],VxAverageEnum); + inputs->GetParameterValue(&vy,&gauss_l1l2l3[0],VyAverageEnum); + + DL_scalar=-dt*gauss_weight*Jdettria; + + DLprime[0][0]=DL_scalar*vx; + DLprime[1][1]=DL_scalar*vy; + + //Do the triple product tL*D*L. + TripleMultiply( &B[0][0],2,numdof,1, + &DLprime[0][0],2,2,0, + &Bprime[0][0],2,numdof,0, + &Ke_gg2[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ + for( i=0; iparameters->FindParam(&dt,DtEnum); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + //recover material parameters + mixed_layer_capacity=matpar->GetMixedLayerCapacity(); + thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + + + GaussTria (&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gauss (nodes on the bedrock) */ + for (ig=0; igGetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); + inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); + + /* Add value into pe_g: */ + for( i=0; iGetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); + inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); + inputs->GetParameterValue(&dhdt_g, &gauss_l1l2l3[0],DhDtEnum); + + /* Add value into pe_g: */ + for( i=0; iGetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); + inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); + + /* Add value into pe_g: */ + for( i=0; iGetParameterValue(&meltingvalue, &gauss_l1l2l3[0],MeltingRateEnum); + + /*Get velocity at gaussian point: */ + inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); + inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); + + /*Get bed slope: */ + inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],BedEnum); + dbdx=slope[0]; + dbdy=slope[1]; + + /* Get Jacobian determinant: */ + GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + //Get L matrix if viscous basal drag present: + GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,NDOF1); + + + /*Build gaussian vector: */ + for(i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&drag_type,DragTypeEnum); + + /*First, if we are on water, return empty vector: */ + if(onwater)return; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + + /* Get gaussian points and weights: */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/ + + /* Start looping on the number of gaussian points: */ + for (ig=0; igGetParameterValue(&thickness, &gauss_l1l2l3[0],ThicknessEnum); + inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],SurfaceEnum); + + /*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the + * element itself: */ + if(drag_type==1){ + inputs->GetParameterValue(&plastic_stress, &gauss_l1l2l3[0],DragCoefficientEnum); + } + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /*Get nodal functions: */ + GetNodalFunctions(l1l2l3, gauss_l1l2l3); + + /*Compute driving stress: */ + driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness; + + /*Build pe_g_gaussian vector: */ + if(drag_type==1){ + for (i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); + + /*If on water, skip: */ + if(onwater)return; + + /*Spawn 3 sing elements: */ + for(i=0;i<3;i++){ + sing=(Sing*)SpawnSing(i); + sing->CreatePVector(pg); + } + + /*clean up*/ + delete sing; + +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorPrognostic {{{1*/ +void Tria::CreatePVectorPrognostic(Vec pg){ + + + /* local declarations */ + int i,j; + + /* node data: */ + const int numgrids=3; + const int NDOF1=1; + const int numdof=NDOF1*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdof]; + int numberofdofspernode; + + /* gaussian points: */ + 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_l1l2l3[3]; + + /* matrix */ + double pe_g[numgrids]={0.0}; + double L[numgrids]; + double Jdettria; + + /*input parameters for structural analysis (diagnostic): */ + double accumulation_g; + double melting_g; + double thickness_g; + + /*parameters: */ + double dt; + + /*retrieve some parameters: */ + this->parameters->FindParam(&dt,DtEnum); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* 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; igGetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); + inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); + inputs->GetParameterValue(&thickness_g, &gauss_l1l2l3[0],ThicknessEnum); + + /* Add value into pe_g: */ + for( i=0; iparameters->FindParam(&dt,DtEnum); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* 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; igGetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); + inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); + inputs->GetParameterValue(&thickness_g, &gauss_l1l2l3[0],ThicknessEnum); + + /* Add value into pe_g: */ + for( i=0; iFindParam(&sub_analysis_type,AnalysisTypeEnum); + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + + /* Get gaussian points and weights: */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/ + + + /* Start looping on the number of gaussian points: */ + for (ig=0; igGetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],SurfaceEnum); + } + if ( (sub_analysis_type==BedSlopeXAnalysisEnum) || (sub_analysis_type==BedSlopeYAnalysisEnum)){ + inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],BedEnum); + } + + /* Get Jacobian determinant: */ + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /*Get nodal functions: */ + GetNodalFunctions(l1l2l3, gauss_l1l2l3); + + /*Build pe_g_gaussian vector: */ + if ( (sub_analysis_type==SurfaceSlopeXAnalysisEnum) || (sub_analysis_type==BedSlopeXAnalysisEnum)){ + for(i=0;iGetMixedLayerCapacity(); + thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + beta=matpar->GetBeta(); + meltingpoint=matpar->GetMeltingPoint(); + + /*retrieve some solution parameters: */ + this->parameters->FindParam(&dt,DtEnum); + + /* Ice/ocean heat exchange flux on ice shelf base */ + + GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gauss 2d (nodes on the bedrock) */ + for (ig=0; igGetParameterValue(&pressure, &gauss_coord[0],PressureEnum); + t_pmp=meltingpoint-beta*pressure; + + /*Calculate scalar parameter*/ + scalar_ocean=gauss_weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice); + if(dt){ + scalar_ocean=dt*scalar_ocean; + } + + for(i=0;i<3;i++){ + P_terms[i]+=scalar_ocean*l1l2l3[i]; + } + } + + /*Add pe_g to global vector pg: */ + VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES); + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + +} +/*}}}*/ +/*FUNCTION Tria::CreatePVectorThermalSheet {{{1*/ +void Tria::CreatePVectorThermalSheet( Vec pg){ + + int i,found; + + const int numgrids=3; + const int NDOF1=1; + const int numdof=numgrids*NDOF1; + int doflist[numdof]; + int numberofdofspernode; + double xyz_list[numgrids][3]; + + double rho_ice; + double heatcapacity; + + /*inputs: */ + double dt; + double pressure_list[3]; + double pressure; + int drag_type; + double basalfriction; + Friction* friction=NULL; + double alpha2,vx,vy; + double geothermalflux_value; + + /* gaussian points: */ + int num_area_gauss,ig; + double* gauss_weights = NULL; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double gauss_weight; + double gauss_coord[3]; + + /*matrices: */ + double Jdet; + double P_terms[numdof]={0.0}; + double l1l2l3[numgrids]; + double scalar; + + int analysis_type; + + /*retrive parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + //recover material parameters + rho_ice=matpar->GetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + + /*retrieve some parameters: */ + this->parameters->FindParam(&dt,DtEnum); + + /*Build frictoin element, needed later: */ + inputs->GetParameterValue(&drag_type,DragTypeEnum); + if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); + friction=new Friction("3d",inputs,matpar,analysis_type); + + /* Ice/ocean heat exchange flux on ice shelf base */ + GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gauss 2d (nodes on the bedrock) */ + for (ig=0; igGetParameterValue(&geothermalflux_value, &gauss_coord[0],GeothermalFluxEnum); + + friction->GetAlpha2(&alpha2,&gauss_coord[0],VxEnum,VyEnum,VzEnum); + inputs->GetParameterValue(&vx, &gauss_coord[0],VxEnum); + inputs->GetParameterValue(&vy, &gauss_coord[0],VyEnum); + basalfriction= alpha2*(pow(vx,(double)2.0)+pow(vy,(double)2.0)); + + /*Calculate scalar parameter*/ + scalar=gauss_weight*Jdet*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice); + if(dt){ + scalar=dt*scalar; + } + + for(i=0;i<3;i++){ + P_terms[i]+=scalar*l1l2l3[i]; + } + } + + /*Add pe_g to global vector pg: */ + VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES); + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + delete friction; + +} +/*}}}*/ +/*FUNCTION Tria::GetArea {{{1*/ +double Tria::GetArea(void){ + + double area=0; + const int numgrids=3; + double xyz_list[numgrids][3]; + double x1,y1,x2,y2,x3,y3; + + /*Get xyz list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + x1=xyz_list[0][0]; y1=xyz_list[0][1]; + x2=xyz_list[1][0]; y2=xyz_list[1][1]; + x3=xyz_list[2][0]; y3=xyz_list[2][1]; + + return x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1; +} +/*}}}*/ +/*FUNCTION Tria::GetAreaCoordinate {{{1*/ +double Tria::GetAreaCoordinate(double x, double y, int which_one){ + + double area=0; + const int numgrids=3; + double xyz_list[numgrids][3]; + double x1,y1,x2,y2,x3,y3; + + /*Get area: */ + area=this->GetArea(); + + /*Get xyz list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + x1=xyz_list[0][0]; y1=xyz_list[0][1]; + x2=xyz_list[1][0]; y2=xyz_list[1][1]; + x3=xyz_list[2][0]; y3=xyz_list[2][1]; + + if(which_one==1){ + /*Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area*/ + return ((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area; + } + else if(which_one==2){ + /*Get second area coordinate = det(x1-x3 x-x3 ; y1-y3 y-y3)/area*/ + return ((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area; + } + else if(which_one==3){ + /*Get third area coordinate 1-area1-area2: */ + return 1-((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area -((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area; + } + else ISSMERROR("%s%i%s\n"," error message: area coordinate ",which_one," done not exist!"); +} +/*}}}*/ +/*FUNCTION Tria::GetB {{{1*/ + +void Tria::GetB(double* B, double* xyz_list, double* gauss_l1l2l3){ + + /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. + * For grid i, Bi can be expressed in the actual coordinate system + * by: + * Bi=[ dh/dx 0 ] + * [ 0 dh/dy ] + * [ 1/2*dh/dy 1/2*dh/dx ] + * where h is the interpolation function for grid i. + * + * We assume B has been allocated already, of size: 3x(NDOF2*numgrids) + */ + + int i; + const int NDOF2=2; + const int numgrids=3; + + double dh1dh3[NDOF2][numgrids]; + + + /*Get dh1dh2dh3 in actual coordinate system: */ + GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss_l1l2l3); + + /*Build B: */ + for (i=0;iGetDofList(&doflist_per_node[0],&numberofdofspernode); + for(j=0;jGetDofList1(); + } + +} +/*}}}*/ +/*FUNCTION Tria::GetJacobian {{{1*/ +void Tria::GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3){ + + /*The Jacobian is constant over the element, discard the gaussian points. + * J is assumed to have been allocated of size NDOF2xNDOF2.*/ + + const int NDOF2=2; const int numgrids=3; - double gaussgrids[numgrids][3]={{1,0,0},{0,1,0},{0,0,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 Tria::MinVx(double* pminvx, bool process_units);{{{1*/ -void Tria::MinVx(double* pminvx, bool process_units){ + + /* node data: */ + const int numgrids=3; + const int NDOF2=2; + double xyz_list[numgrids][3]; + int doflist1[numgrids]; + double dh1dh3[NDOF2][numgrids]; + + /* grid data: */ + double drag; + double alpha_complement; + Friction* friction=NULL; + + /* gaussian points: */ + 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_l1l2l3[3]; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; + + /* parameters: */ + double vx,vy,vz; + double lambda,mu,xi; + double bed,thickness,Neff; + double surface_normal[3]; + double bed_normal[3]; + double dk[NDOF2]; + + /*element vector at the gaussian points: */ + double grade_g[numgrids]={0.0}; + double grade_g_gaussian[numgrids]; + + /* Jacobian: */ + double Jdet; + + /*nodal functions: */ + double l1l2l3[3]; + + /* strain rate: */ + double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + + /*inputs: */ + bool shelf; + int drag_type; + + /*parameters: */ + double cm_noisedmp; + double cm_mindmp_slope; + double cm_mindmp_value; + double cm_maxdmp_value; + double cm_maxdmp_slope; + + int analysis_type; + + /*retrive parameters: */ + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + + /*retrieve inputs :*/ + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + inputs->GetParameterValue(&drag_type,DragTypeEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); + this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum); + this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum); + this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum); + this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum); + + /*Get out if shelf*/ + if(shelf)return; + + /* Get node coordinates and dof list: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList1(&doflist1[0]); + + /*Build frictoin element, needed later: */ + inputs->GetParameterValue(&drag_type,DragTypeEnum); + friction=new Friction("2d",inputs,matpar,analysis_type); + + /* 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, 4); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igGetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum); + else alpha_complement=0; + inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum); + + /*recover lambda mu and xi: */ + inputs->GetParameterValue(&lambda, &gauss_l1l2l3[0],AdjointxEnum); + inputs->GetParameterValue(&mu, &gauss_l1l2l3[0],AdjointyEnum); + inputs->GetParameterValue(&xi, &gauss_l1l2l3[0],AdjointzEnum); + + /*recover vx vy and vz: */ + inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); + inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); + inputs->GetParameterValue(&vz, &gauss_l1l2l3[0],VzEnum); + + /*Get normal vector to the bed */ + SurfaceNormal(&surface_normal[0],xyz_list); + + 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]; + + /* Get Jacobian determinant: */ + GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /* Get nodal functions value at gaussian point:*/ + GetNodalFunctions(l1l2l3, gauss_l1l2l3); + + /*Get nodal functions derivatives*/ + GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3); + + /*Get k derivative: dk/dx */ + inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); + + /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ + for (i=0;icm_maxdmp_value){ + grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; + } + } + + /*Add gradje_g_gaussian vector to gradje_g: */ + for( i=0; iinputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0])); + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + delete friction; + +} +/*}}}*/ +/*FUNCTION Tria::SetClone {{{1*/ +void Tria::SetClone(int* minranks){ + + ISSMERROR("not implemented yet"); +} +/*}}}1*/ +/*FUNCTION Tria::SpawnBeam {{{1*/ +void* Tria::SpawnBeam(int g0, int g1){ int i; - int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,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;iparameters; + beam_inputs=(Inputs*)this->inputs->SpawnBeamInputs(indices); + + beam=new Beam(); + beam->id=this->id; + beam->inputs=beam_inputs; + beam->parameters=beam_parameters; + + /*now deal with nodes, matice and matpar: */ + beam->nodes=(Node**)xmalloc(2*sizeof(Node*)); + for(i=0;i<2;i++)beam->nodes[i]=this->nodes[indices[i]]; + beam->matice=this->matice; + beam->matpar=this->matpar; + + + return beam; +} +/*}}}*/ +/*FUNCTION Tria::SpawnSing {{{1*/ +void* Tria::SpawnSing(int index){ + + Sing* sing=NULL; + int zero=0; + Parameters *sing_parameters = NULL; + Inputs *sing_inputs = NULL; + + sing_parameters=this->parameters; + sing_inputs=(Inputs*)this->inputs->SpawnSingInputs(index); + + sing=new Sing(); + sing->id=this->id; + sing->inputs=sing_inputs; + sing->parameters=sing_parameters; + + /*now deal with node,matice and matpar: */ + sing->node=this->nodes[index]; + sing->matice=this->matice; + sing->matpar=this->matpar; + + return sing; +} +/*}}}*/ +/*FUNCTION Tria::SurfaceNormal{{{1*/ + +void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){ int i; - int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,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 Tria::MaxAbsVx(double* pmaxabsvx, bool process_units);{{{1*/ -void Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){ - - int i; - int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,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 Tria::MinVy(double* pminvy, bool process_units);{{{1*/ -void Tria::MinVy(double* pminvy, bool process_units){ - - int i; - int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,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 Tria::MaxAbsVy(double* pmaxabsvy, bool process_units);{{{1*/ -void Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){ - - int i; - int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,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 Tria::MinVz(double* pminvz, bool process_units);{{{1*/ -void Tria::MinVz(double* pminvz, bool process_units){ - - int i; - int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,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 Tria::MaxAbsVz(double* pmaxabsvz, bool process_units);{{{1*/ -void Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){ - - int i; - int dim; - const int numgrids=3; - double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,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 Tria::InputDuplicate(int original_enum,int new_enum){{{1*/ -void Tria::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 Tria::InputScale(int enum_type,double scale_factor){{{1*/ -void Tria::InputScale(int enum_type,double scale_factor){ - - Input* input=NULL; - - /*Make a copy of the original input: */ - input=(Input*)this->inputs->GetInput(enum_type); - - /*Scale: */ - input->Scale(scale_factor); -} -/*}}}*/ -/*FUNCTION Tria::InputAXPY(int YEnum, double scalar, int XEnum);{{{1*/ -void Tria::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!"); - - /*Scale: */ - yinput->AXPY(xinput,scalar); -} -/*}}}*/ -/*FUNCTION Tria::InputControlConstrain(int control_type, double cm_min, double cm_max){{{1*/ -void Tria::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 Tria::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/ -void Tria::GetVectorFromInputs(Vec vector,int NameEnum){ - - int i; - const int numvertices=3; - 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 Tria::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/ -void Tria::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; - -} -/*}}}*/ + double v13[3]; + double v23[3]; + double normal[3]; + double normal_norm; + + for (i=0;i<3;i++){ + v13[i]=xyz_list[0][i]-xyz_list[2][i]; + v23[i]=xyz_list[1][i]-xyz_list[2][i]; + } + + normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; + normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; + normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; + + normal_norm=sqrt( pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2) ); + + *(surface_normal)=normal[0]/normal_norm; + *(surface_normal+1)=normal[1]/normal_norm; + *(surface_normal+2)=normal[2]/normal_norm; + +} +/*}}}*/ +/*FUNCTION Tria::UpdateFromDakota {{{1*/ +void Tria::UpdateFromDakota(void* vinputs){ + + int i; + int dofs[1]={0}; + double temperature_list[3]; + double temperature_average; + double B_list[3]; + double B_average; + double new_h[3]; + + /*Update internal data if inputs holds new values: */ + /*inputs->Recover("thickness",&this->properties.h[0],1,dofs,3,(void**)nodes); + if(inputs->Recover("thickness",&new_h[0],1,dofs,3,(void**)nodes)){ + //density, needed later: + double di=(this->matpar->GetRhoIce()/this->matpar->GetRhoWater()); + //Go through grids: + for (i=0;i<3;i++){ + if(nodes[i]->IsOnShelf()){ + this->b[i]=this->b[i]-di*(new_h[i]-h[i]); //hydrostatic equilibrium; + } + this->s[i]=this->b[i]+new_h[i]; + this->h[i]=new_h[i]; + } + }*/ + + ISSMERROR("not supported yet!"); + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 4284) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 4285) @@ -54,11 +54,4 @@ /*Update virtual functions resolution: {{{1*/ void InputUpdateFromSolution(double* solutiong); - void InputUpdateFromSolutionDiagnosticHoriz( double* solution); - void InputUpdateFromSolutionSlopeCompute( double* solution); - void InputUpdateFromSolutionPrognostic( double* solution); - void InputUpdateFromSolutionPrognostic2(double* solution); - void InputUpdateFromSolutionBalancedthickness( double* solution); - void InputUpdateFromSolutionBalancedthickness2( double* solution); - void InputUpdateFromSolutionBalancedvelocities( double* solution); void InputUpdateFromVector(double* vector, int name, int type); void InputUpdateFromVector(int* vector, int name, int type); @@ -67,5 +60,4 @@ void InputUpdateFromConstant(int constant, int name); void InputUpdateFromConstant(bool constant, int name); - void UpdateFromDakota(void* inputs); /*}}}*/ /*Element virtual functions definitions: {{{1*/ @@ -160,4 +152,11 @@ void GetSolutionFromInputsDiagnosticHoriz(Vec solution); void GradjDragStokes(Vec gradient); + void InputUpdateFromSolutionDiagnosticHoriz( double* solution); + void InputUpdateFromSolutionSlopeCompute( double* solution); + void InputUpdateFromSolutionPrognostic( double* solution); + void InputUpdateFromSolutionPrognostic2(double* solution); + void InputUpdateFromSolutionBalancedthickness( double* solution); + void InputUpdateFromSolutionBalancedthickness2( double* solution); + void InputUpdateFromSolutionBalancedvelocities( double* solution); bool IsInput(int name); void SetClone(int* minranks); @@ -165,4 +164,5 @@ void* SpawnSing(int g0); void SurfaceNormal(double* surface_normal, double xyz_list[3][3]); + void UpdateFromDakota(void* inputs); /*}}}*/ Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 4284) +++ /issm/trunk/src/c/objects/Node.cpp (revision 4285) @@ -151,5 +151,6 @@ analysis_type==PrognosticAnalysisEnum || analysis_type==MeltingAnalysisEnum || - analysis_type==SlopeAnalysisEnum || + analysis_type==BedSlopeAnalysisEnum || + analysis_type==SurfaceSlopeAnalysisEnum || analysis_type==BalancedvelocitiesAnalysisEnum || analysis_type==BalancedthicknessAnalysisEnum Index: /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp =================================================================== --- /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp (revision 4284) +++ /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp (revision 4285) @@ -25,5 +25,5 @@ numdofs=2; } - else if (analysis_type==SlopeAnalysisEnum){ + else if (analysis_type==BedSlopeAnalysisEnum || analysis_type==SurfaceSlopeAnalysisEnum){ numdofs=1; } Index: /issm/trunk/src/c/solutions/adjoint_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/adjoint_core.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/adjoint_core.cpp (revision 4285) @@ -20,5 +20,7 @@ char* solverstring=NULL; bool isstokes=false; - bool conserve_loads=true; + bool conserve_loads=true; + int dim; + int solution_type; /*intermediary: */ @@ -37,4 +39,6 @@ femmodel->parameters->FindParam(&verbose,VerboseEnum); femmodel->parameters->FindParam(&isstokes,IsStokesEnum); + femmodel->parameters->FindParam(&dim,DimEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); /*set analysis type to compute velocity: */ @@ -72,4 +76,11 @@ InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,adjoint_g); + if(verbose)_printf_("saving results:\n"); + if(solution_type==AdjointAnalysisEnum){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointxEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointyEnum); + if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointzEnum); + } + /*Free ressources:*/ xfree((void**)&solverstring); Index: /issm/trunk/src/c/solutions/balancedthickness2_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedthickness2_core.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/balancedthickness2_core.cpp (revision 4285) @@ -17,4 +17,5 @@ int verbose=0; int dim; + int solution_type; /*activate formulation: */ @@ -24,4 +25,5 @@ femmodel->parameters->FindParam(&verbose,VerboseEnum); femmodel->parameters->FindParam(&dim,DimEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); _printf_("depth averaging velocity...\n"); @@ -42,5 +44,8 @@ if(verbose)_printf_("saving results:\n"); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + if(solution_type==Balancedthickness2AnalysisEnum){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + } + } Index: /issm/trunk/src/c/solutions/balancedthickness_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedthickness_core.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/balancedthickness_core.cpp (revision 4285) @@ -17,4 +17,5 @@ int verbose=0; int dim; + int solution_type; /*activate formulation: */ @@ -24,4 +25,5 @@ femmodel->parameters->FindParam(&verbose,VerboseEnum); femmodel->parameters->FindParam(&dim,DimEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); _printf_("depth averaging velocity...\n"); @@ -37,5 +39,7 @@ if(verbose)_printf_("saving results:\n"); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + if(solution_type==BalancedthicknessAnalysisEnum){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + } } Index: /issm/trunk/src/c/solutions/balancedvelocities_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/balancedvelocities_core.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/balancedvelocities_core.cpp (revision 4285) @@ -16,4 +16,5 @@ int verbose=0; int dim; + int solution_type; /*activate formulation: */ @@ -23,4 +24,5 @@ femmodel->parameters->FindParam(&verbose,VerboseEnum); femmodel->parameters->FindParam(&dim,DimEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); _printf_("depth averaging velocity...\n"); @@ -37,6 +39,8 @@ if(verbose)_printf_("saving results:\n"); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); + if(solution_type==BalancedvelocitiesAnalysisEnum){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); + } } Index: /issm/trunk/src/c/solutions/bedslope.cpp =================================================================== --- /issm/trunk/src/c/solutions/bedslope.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/bedslope.cpp (revision 4285) @@ -37,6 +37,6 @@ double start_init, finish_init; - int analyses[1]={SlopeAnalysisEnum}; - int solution_type=SlopeAnalysisEnum; + int analyses[1]={BedSlopeAnalysisEnum}; + int solution_type=BedSlopeAnalysisEnum; MODULEBOOT(); Index: /issm/trunk/src/c/solutions/bedslope_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/bedslope_core.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/bedslope_core.cpp (revision 4285) @@ -13,12 +13,14 @@ /*parameters: */ - int verbose; - int dim; + int verbose; + int dim; bool isstokes; bool ishutter; + int solution_type; /*Recover some parameters: */ femmodel->parameters->FindParam(&verbose,VerboseEnum); femmodel->parameters->FindParam(&dim,DimEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); if(verbose)_printf_("%s\n","computing slope..."); @@ -38,6 +40,8 @@ if(verbose)_printf_("saving results:\n"); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeXEnum); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeYEnum); + if(solution_type==BedSlopeAnalysisEnum){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeXEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeYEnum); + } } Index: /issm/trunk/src/c/solutions/control_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/control_core.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/control_core.cpp (revision 4285) @@ -26,4 +26,5 @@ double cm_max; int cm_gradient; + int dim; double* fit=NULL; @@ -57,4 +58,5 @@ femmodel->parameters->FindParam(&cm_gradient,CmGradientEnum); femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum); + femmodel->parameters->FindParam(&dim,DimEnum); /*}}}*/ @@ -116,4 +118,11 @@ /*some results not computed by steadystate_core or diagnostic_core: */ InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself! + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VelEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,GradientEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointxEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointyEnum); + if(dim==3)InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum); femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,0)); femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumAsString(control_type),1,0)); Index: /issm/trunk/src/c/solutions/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 4285) @@ -24,4 +24,5 @@ bool conserve_loads=true; bool modify_loads=true; + int solution_type; /* recover parameters:*/ @@ -33,4 +34,5 @@ femmodel->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); /*for qmu analysis, reinitialize velocity so that fake sensitivities do not show up as a result of a different restart of the convergence at each trial.*/ @@ -85,9 +87,10 @@ if(verbose)_printf_("saving results:\n"); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VelEnum); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum); - if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum); - + if(solution_type==DiagnosticAnalysisEnum){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VelEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum); + if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum); + } } Index: /issm/trunk/src/c/solutions/prognostic2_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/prognostic2_core.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/prognostic2_core.cpp (revision 4285) @@ -15,4 +15,5 @@ /*flags: */ int verbose=0; + int solution_type; /*activate formulation: */ @@ -21,4 +22,5 @@ /*recover parameters: */ femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); _printf_("depth averaging velocity...\n"); @@ -35,5 +37,7 @@ if(verbose)_printf_("saving results:\n"); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + if(solution_type==Prognostic2AnalysisEnum){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + } } Index: /issm/trunk/src/c/solutions/prognostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/prognostic_core.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/prognostic_core.cpp (revision 4285) @@ -15,4 +15,5 @@ /*parameters: */ int verbose=0; + int solution_type; /*activate formulation: */ @@ -21,4 +22,5 @@ /*recover parameters: */ femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); _printf_("depth averaging velocity...\n"); @@ -33,5 +35,7 @@ if(verbose)_printf_("saving results:\n"); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + if(solution_type==PrognosticAnalysisEnum){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + } } Index: /issm/trunk/src/c/solutions/steadystate_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/steadystate_core.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/steadystate_core.cpp (revision 4285) @@ -20,8 +20,10 @@ int verbose; int dim; + int solution_type; /* recover parameters:*/ femmodel->parameters->FindParam(&verbose,VerboseEnum); femmodel->parameters->FindParam(&dim,DimEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); /*intialize counters: */ @@ -56,8 +58,10 @@ if(verbose)_printf_("saving results:\n"); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum); - if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum); + if(solution_type==SteadystateAnalysisEnum){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum); + if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum); + } } Index: /issm/trunk/src/c/solutions/surfaceslope.cpp =================================================================== --- /issm/trunk/src/c/solutions/surfaceslope.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/surfaceslope.cpp (revision 4285) @@ -37,6 +37,6 @@ double start_init, finish_init; - int analyses[1]={SlopeAnalysisEnum}; - int solution_type=SlopeAnalysisEnum; + int analyses[1]={SurfaceSlopeAnalysisEnum}; + int solution_type=SurfaceSlopeAnalysisEnum; MODULEBOOT(); Index: /issm/trunk/src/c/solutions/surfaceslope_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/surfaceslope_core.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/surfaceslope_core.cpp (revision 4285) @@ -17,8 +17,10 @@ bool isstokes; bool ishutter; + int solution_type; /*Recover some parameters: */ femmodel->parameters->FindParam(&verbose,VerboseEnum); femmodel->parameters->FindParam(&dim,DimEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); if(verbose)_printf_("%s\n","computing slope..."); @@ -38,6 +40,8 @@ if(verbose)_printf_("saving results:\n"); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeXEnum); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeYEnum); + if(solution_type==SurfaceSlopeAnalysisEnum){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeXEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeYEnum); + } } Index: /issm/trunk/src/c/solutions/thermal_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/thermal_core.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/thermal_core.cpp (revision 4285) @@ -23,4 +23,5 @@ double dt; double melting_offset; + int solution_type; //first recover parameters common to all solutions @@ -28,4 +29,5 @@ femmodel->parameters->FindParam(&ndt,NdtEnum); femmodel->parameters->FindParam(&dt,DtEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); /*Compute number of time steps: */ @@ -46,6 +48,8 @@ if(verbose)_printf_("saving results:\n"); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,i,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingRateEnum,i,time); + if(solution_type==ThermalAnalysisEnum){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,i,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingRateEnum,i,time); + } } Index: /issm/trunk/src/c/solutions/transient2d_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient2d_core.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/transient2d_core.cpp (revision 4285) @@ -18,4 +18,5 @@ double finaltime; double dt; + int solution_type; /*intermediary: */ @@ -27,4 +28,5 @@ femmodel->parameters->FindParam(&finaltime,NdtEnum); femmodel->parameters->FindParam(&dt,DtEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); /*initialize: */ @@ -36,6 +38,7 @@ _printf_("%s%g%s%i%s%g\n","time [yr]: ",time," iteration number: ",step,"/",floor(finaltime/dt)); + /*Set step and time: */ + time+=dt; step+=1; - time+=dt; if(verbose)_printf_("%s\n","computing new velocity"); @@ -49,10 +52,13 @@ if(verbose)_printf_("%s\n","saving results:\n"); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); + if(solution_type==Transient2DAnalysisEnum){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VelEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); + } } Index: /issm/trunk/src/c/solutions/transient3d_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/transient3d_core.cpp (revision 4284) +++ /issm/trunk/src/c/solutions/transient3d_core.cpp (revision 4285) @@ -18,4 +18,5 @@ double finaltime; double dt; + int solution_type; /*intermediary: */ @@ -27,4 +28,5 @@ femmodel->parameters->FindParam(&finaltime,NdtEnum); femmodel->parameters->FindParam(&dt,DtEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); /*initialize: */ @@ -58,13 +60,15 @@ if(verbose)_printf_("%s\n","saving results:\n"); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time); - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingRateEnum,step,time); + if(solution_type==Transient3DAnalysisEnum){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingRateEnum,step,time); + } if (step%5==0){