Ignore:
Timestamp:
06/28/10 23:22:14 (15 years ago)
Author:
Eric.Larour
Message:

Cleanup of Tria, Sing and Beam elements + New flags to stop cores from creating too many results

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/Sing.cpp

    r4248 r4285  
    3131}
    3232/*}}}*/
    33 /*FUNCTION void Sing::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);{{{1*/
    34 void Sing::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){
    35         ISSMERROR(" not supported yet!");
    36 }
    37 /*}}}*/
    3833
    3934/*Object virtual functions definitions: */
     35/*FUNCTION Sing::copy {{{1*/
     36Object* Sing::copy() {
     37
     38        int i;
     39        Sing* sing=NULL;
     40
     41        sing=new Sing();
     42
     43        /*copy fields: */
     44        sing->id=this->id;
     45        if(this->inputs){
     46                sing->inputs=(Inputs*)this->inputs->Copy();
     47        }
     48        else{
     49                sing->inputs=new Inputs();
     50        }
     51        /*point parameters: */
     52        sing->parameters=this->parameters;
     53
     54        /*pointers: */
     55        sing->node=this->node;
     56        sing->matice=this->matice;
     57        sing->matpar=this->matpar;
     58
     59        return sing;
     60}
     61/*}}}*/
     62/*FUNCTION Sing::DeepEcho {{{1*/
     63void Sing::DeepEcho(void){
     64
     65        printf("Sing:\n");
     66        printf("   id: %i\n",id);
     67        node->DeepEcho();
     68        matice->DeepEcho();
     69        matpar->DeepEcho();
     70        printf("   parameters\n");
     71        parameters->DeepEcho();
     72        printf("   inputs\n");
     73        inputs->DeepEcho();
     74
     75        return;
     76}
     77/*}}}*/
     78/*FUNCTION Sing::Demarshall {{{1*/
     79void  Sing::Demarshall(char** pmarshalled_dataset){
     80        ISSMERROR(" not supported yet!");
     81}
     82/*}}}*/
    4083/*FUNCTION Sing::Echo{{{1*/
    4184
    … …  
    5396}
    5497/*}}}*/
    55 /*FUNCTION Sing::DeepEcho {{{1*/
    56 void Sing::DeepEcho(void){
    57 
    58         printf("Sing:\n");
    59         printf("   id: %i\n",id);
    60         node->DeepEcho();
    61         matice->DeepEcho();
    62         matpar->DeepEcho();
    63         printf("   parameters\n");
    64         parameters->DeepEcho();
    65         printf("   inputs\n");
    66         inputs->DeepEcho();
    67 
    68         return;
     98/*FUNCTION Sing::Enum {{{1*/
     99int Sing::Enum(void){
     100
     101        return SingEnum;
     102
    69103}
    70104/*}}}*/
    71105/*FUNCTION Sing::Id {{{1*/
    72106int    Sing::Id(void){ return id; }
     107/*}}}*/
     108/*FUNCTION Sing::Marshall {{{1*/
     109void  Sing::Marshall(char** pmarshalled_dataset){
     110        ISSMERROR(" not supported yet!");
     111}
     112/*}}}*/
     113/*FUNCTION Sing::MashallSize {{{1*/
     114int   Sing::MarshallSize(){
     115        ISSMERROR(" not supported yet!");
     116}
    73117/*}}}*/
    74118/*FUNCTION Sing::MyRank {{{1*/
    … …  
    78122}
    79123/*}}}*/
    80 /*FUNCTION Sing::Marshall {{{1*/
    81 void  Sing::Marshall(char** pmarshalled_dataset){
    82         ISSMERROR(" not supported yet!");
    83 }
    84 /*}}}*/
    85 /*FUNCTION Sing::MashallSize {{{1*/
    86 int   Sing::MarshallSize(){
    87         ISSMERROR(" not supported yet!");
    88 }
    89 /*}}}*/
    90 /*FUNCTION Sing::Demarshall {{{1*/
    91 void  Sing::Demarshall(char** pmarshalled_dataset){
    92         ISSMERROR(" not supported yet!");
    93 }
    94 /*}}}*/
    95 /*FUNCTION Sing::Enum {{{1*/
    96 int Sing::Enum(void){
    97 
    98         return SingEnum;
    99 
    100 }
    101 /*}}}*/
    102 /*FUNCTION Sing::copy {{{1*/
    103 Object* Sing::copy() {
    104 
    105         int i;
    106         Sing* sing=NULL;
    107 
    108         sing=new Sing();
    109 
    110         /*copy fields: */
    111         sing->id=this->id;
    112         if(this->inputs){
    113                 sing->inputs=(Inputs*)this->inputs->Copy();
    114         }
    115         else{
    116                 sing->inputs=new Inputs();
    117         }
    118         /*point parameters: */
    119         sing->parameters=this->parameters;
    120 
    121         /*pointers: */
    122         sing->node=this->node;
    123         sing->matice=this->matice;
    124         sing->matpar=this->matpar;
    125 
    126         return sing;
    127 }
    128 /*}}}*/
    129 
    130 /*Sing management*/
     124
     125/*Update virtual functions definitions: */
     126/*FUNCTION Sing::InputUpdateFromSolution {{{1*/
     127void  Sing::InputUpdateFromSolution(double* solution){
     128        ISSMERROR(" not supported yet!");
     129}
     130/*}}}*/
     131/*FUNCTION Sing::InputUpdateFromVector(double* vector, int name, int type);{{{1*/
     132void  Sing::InputUpdateFromVector(double* vector, int name, int type){
     133
     134        /*Check that name is an element input*/
     135        if (!IsInput(name)) return;
     136
     137        switch(type){
     138
     139                case VertexEnum:
     140
     141                        /*New SingVertexInpu*/
     142                        double value;
     143
     144                        /*Get values on the 6 vertices*/
     145                        value=vector[node->GetVertexDof()];
     146
     147                        /*update input*/
     148                        this->inputs->AddInput(new SingVertexInput(name,value));
     149                        return;
     150
     151                default:
     152                        ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type));
     153        }
     154}
     155/*}}}*/
     156/*FUNCTION Sing::InputUpdateFromVector(int* vector, int name, int type);{{{1*/
     157void  Sing::InputUpdateFromVector(int* vector, int name, int type){
     158        ISSMERROR(" not supported yet!");
     159}
     160/*}}}*/
     161/*FUNCTION Sing::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/
     162void  Sing::InputUpdateFromVector(bool* vector, int name, int type){
     163        ISSMERROR(" not supported yet!");
     164}
     165/*}}}*/
     166
     167/*Element virtual functions definitions: */
     168/*FUNCTION Sing::ComputeBasalStress {{{1*/
     169void  Sing::ComputeBasalStress(Vec p_g){
     170
     171        ISSMERROR("Not implemented yet");
     172
     173}
     174/*}}}*/
     175/*FUNCTION Sing::ComputePressure {{{1*/
     176void  Sing::ComputePressure(Vec p_g){
     177
     178        int    dof;
     179        double pressure;
     180        double thickness;
     181        double rho_ice,g;
     182
     183        /*Get dof list on which we will plug the pressure values: */
     184        GetDofList1(&dof);
     185
     186        /*pressure is lithostatic: */
     187        rho_ice=matpar->GetRhoIce();
     188        g=matpar->GetG();
     189        inputs->GetParameterValue(&thickness,ThicknessEnum);
     190        pressure=rho_ice*g*thickness;
     191       
     192        /*plug local pressure values into global pressure vector: */
     193        VecSetValue(p_g,dof,pressure,INSERT_VALUES);
     194
     195}
     196/*}}}*/
     197/*FUNCTION Sing::ComputeStrainRate {{{1*/
     198void  Sing::ComputeStrainRate(Vec p_g){
     199
     200        ISSMERROR("Not implemented yet");
     201
     202}
     203/*}}}*/
    131204/*FUNCTION Sing::Configure {{{1*/
    132205void  Sing::Configure(Elements* elementsin,Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){
    … …  
    136209}
    137210/*}}}*/
     211/*FUNCTION Sing::CostFunction {{{1*/
     212double Sing::CostFunction(){
     213        ISSMERROR(" not supported yet!");
     214}
     215/*}}}*/
     216/*FUNCTION Sing::CreateKMatrix {{{1*/
     217
     218void  Sing::CreateKMatrix(Mat Kgg){
     219
     220        int analysis_type;
     221
     222        /*retrive parameters: */
     223        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     224
     225        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
     226        if (analysis_type==DiagnosticHutterAnalysisEnum){
     227                CreateKMatrixDiagnosticHutter( Kgg);
     228        }
     229        else{
     230                ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
     231        }
     232
     233}
     234/*}}}*/
     235/*FUNCTION Sing::CreatePVector {{{1*/
     236void  Sing::CreatePVector(Vec pg){
     237
     238        int analysis_type;
     239
     240        /*retrive parameters: */
     241        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     242       
     243        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
     244        if (analysis_type==DiagnosticHutterAnalysisEnum){
     245                        CreatePVectorDiagnosticHutter( pg);
     246        }
     247        else{
     248                ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
     249        }
     250
     251}
     252/*}}}*/
     253/*FUNCTION Sing::Du {{{1*/
     254void  Sing::Du(Vec){
     255        ISSMERROR(" not supported yet!");
     256}
     257/*}}}*/
     258/*FUNCTION Sing::GetBedList {{{1*/
     259void  Sing::GetBedList(double*){
     260        ISSMERROR(" not supported yet!");
     261}
     262/*}}}*/
     263/*FUNCTION Sing::GetMatPar {{{1*/
     264void* Sing::GetMatPar(){
     265
     266        return matpar;
     267}
     268/*}}}*/
     269/*FUNCTION Sing::GetNodes {{{1*/
     270void  Sing::GetNodes(void** vpnodes){
     271       
     272        Node** pnodes=NULL;
     273
     274        /*recover nodes: */
     275        pnodes=(Node**)vpnodes;
     276
     277        pnodes[0]=node;
     278}
     279/*}}}*/
     280/*FUNCTION Sing::GetOnBed {{{1*/
     281bool   Sing::GetOnBed(){
     282        ISSMERROR(" not supported yet!");
     283}
     284/*}}}*/
     285/*FUNCTION Sing::GetShelf {{{1*/
     286bool   Sing::GetShelf(){
     287        ISSMERROR(" not supported yet!");
     288}
     289/*}}}*/
     290/*FUNCTION Sing::GetSolutionFromInputs(Vec solution);{{{1*/
     291void  Sing::GetSolutionFromInputs(Vec solution){
     292        ISSMERROR(" not supported yet!");
     293}
     294/*}}}*/
     295/*FUNCTION Sing::GetThicknessList {{{1*/
     296void  Sing::GetThicknessList(double* thickness_list){
     297        ISSMERROR(" not supported yet!");
     298}
     299/*}}}*/
     300/*FUNCTION Sing::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/
     301void  Sing::GetVectorFromInputs(Vec vector,int NameEnum){
     302
     303        int i;
     304        const int numvertices=1;
     305        int doflist1[numvertices];
     306
     307        /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
     308        for(i=0;i<this->inputs->Size();i++){
     309                Input* input=(Input*)this->inputs->GetObjectByOffset(i);
     310                if(input->EnumType()==NameEnum){
     311                        /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
     312                        this->GetDofList1(&doflist1[0]);
     313                        input->GetVectorFromInputs(vector,&doflist1[0]);
     314                        break;
     315                }
     316        }
     317}
     318/*}}}*/
     319/*FUNCTION Sing::Gradj {{{1*/
     320void  Sing::Gradj(Vec gradient,int control_type){
     321        ISSMERROR(" not supported yet!");
     322}
     323/*}}}*/
     324/*FUNCTION Sing::GradB {{{1*/
     325void  Sing::GradjB(Vec gradient){
     326        ISSMERROR(" not supported yet!");
     327}
     328/*}}}*/
     329/*FUNCTION Sing::GradjDrag {{{1*/
     330void  Sing::GradjDrag(Vec gradient){
     331        ISSMERROR(" not supported yet!");
     332}
     333/*}}}*/
     334/*FUNCTION Sing::InputAXPY(int YEnum, double scalar, int XEnum);{{{1*/
     335void  Sing::InputAXPY(int YEnum, double scalar, int XEnum){
     336
     337        Input* xinput=NULL;
     338        Input* yinput=NULL;
     339
     340        /*Find x and y inputs: */
     341        xinput=(Input*)this->inputs->GetInput(XEnum);
     342        yinput=(Input*)this->inputs->GetInput(YEnum);
     343
     344        /*some checks: */
     345        if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!");
     346        if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!");
     347
     348        /*Scale: */
     349        yinput->AXPY(xinput,scalar);
     350}
     351/*}}}*/
     352/*FUNCTION Sing::InputControlConstrain(int control_type, double cm_min, double cm_max){{{1*/
     353void  Sing::InputControlConstrain(int control_type, double cm_min, double cm_max){
     354
     355        Input* input=NULL;
     356
     357        /*Find input: */
     358        input=(Input*)this->inputs->GetInput(control_type);
     359       
     360        /*Do nothing if we  don't find it: */
     361        if(!input)return;
     362
     363        /*Constrain input using cm_min and cm_max: */
     364        input->Constrain(cm_min,cm_max);
     365
     366}
     367/*}}}*/
     368/*FUNCTION Sing::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/
     369void  Sing::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){
     370
     371        int i;
     372        Input** new_inputs=NULL;
     373        Input** old_inputs=NULL;
     374        int     converged=1;
     375
     376        new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
     377        old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
     378       
     379        for(i=0;i<num_enums/2;i++){
     380                new_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+0]);
     381                old_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+1]);
     382                if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(enums[2*i+0]));
     383                if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(enums[2*i+0]));
     384        }
     385
     386        /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
     387        for(i=0;i<num_criterionenums;i++){
     388                IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,criterionenums[i]);
     389                if(eps[i]>criterionvalues[i]) converged=0;
     390        }
     391
     392        /*Assign output pointers:*/
     393        *pconverged=converged;
     394
     395}
     396/*}}}*/
     397/*FUNCTION Sing::InputDuplicate(int original_enum,int new_enum){{{1*/
     398void  Sing::InputDuplicate(int original_enum,int new_enum){
     399
     400        Input* original=NULL;
     401        Input* copy=NULL;
     402
     403        /*Make a copy of the original input: */
     404        original=(Input*)this->inputs->GetInput(original_enum);
     405        copy=(Input*)original->copy();
     406
     407        /*Change copy enum to reinitialized_enum: */
     408        copy->ChangeEnum(new_enum);
     409
     410        /*Add copy into inputs, it will wipe off the one already there: */
     411        inputs->AddObject((Input*)copy);
     412}
     413/*}}}*/
     414/*FUNCTION Sing::InputScale(int enum_type,double scale_factor){{{1*/
     415void  Sing::InputScale(int enum_type,double scale_factor){
     416
     417        Input* input=NULL;
     418
     419        /*Make a copy of the original input: */
     420        input=(Input*)this->inputs->GetInput(enum_type);
     421
     422        /*Scale: */
     423        input->Scale(scale_factor);
     424}
     425/*}}}*/
     426/*FUNCTION Sing::InputToResult(int enum_type,int step,double time){{{1*/
     427void  Sing::InputToResult(int enum_type,int step,double time){
     428        ISSMERROR(" not supported yet!");
     429}
     430/*}}}*/
     431/*FUNCTION Sing::MassFlux {{{1*/
     432double Sing::MassFlux( double* segment){
     433        ISSMERROR(" not supported yet!");
     434}
     435/*}}}*/
     436/*FUNCTION Sing::MaxAbsVx(double* pmaxabsvx, bool process_units);{{{1*/
     437void  Sing::MaxAbsVx(double* pmaxabsvx, bool process_units){
     438
     439        int dim;
     440        double  maxabsvx;
     441
     442        /*retrieve dim parameter: */
     443        parameters->FindParam(&dim,DimEnum);
     444
     445        /*retrive velocity values at nodes */
     446        inputs->GetParameterValue(&maxabsvx,VxEnum);
     447        maxabsvx=fabs(maxabsvx);
     448
     449        /*Assign output pointers:*/
     450        *pmaxabsvx=maxabsvx;
     451}
     452/*}}}*/
     453/*FUNCTION Sing::MaxAbsVy(double* pmaxabsvy, bool process_units);{{{1*/
     454void  Sing::MaxAbsVy(double* pmaxabsvy, bool process_units){
     455
     456        int dim;
     457        double  maxabsvy;
     458
     459        /*retrieve dim parameter: */
     460        parameters->FindParam(&dim,DimEnum);
     461
     462        /*retrive velocity values at nodes */
     463        inputs->GetParameterValue(&maxabsvy,VyEnum);
     464        maxabsvy=fabs(maxabsvy);
     465
     466        /*Assign output pointers:*/
     467        *pmaxabsvy=maxabsvy;
     468}
     469/*}}}*/
     470/*FUNCTION Sing::MaxAbsVz(double* pmaxabsvz, bool process_units);{{{1*/
     471void  Sing::MaxAbsVz(double* pmaxabsvz, bool process_units){
     472
     473        int dim;
     474        double  maxabsvz;
     475
     476        /*retrieve dim parameter: */
     477        parameters->FindParam(&dim,DimEnum);
     478
     479        /*retrive velocity values at nodes */
     480        inputs->GetParameterValue(&maxabsvz,VzEnum);
     481        maxabsvz=fabs(maxabsvz);
     482
     483        /*Assign output pointers:*/
     484        *pmaxabsvz=maxabsvz;
     485}
     486/*}}}*/
     487/*FUNCTION Sing::MaxVel(double* pmaxvel, bool process_units);{{{1*/
     488void  Sing::MaxVel(double* pmaxvel, bool process_units){
     489
     490        int dim;
     491        double  vx;
     492        double  vy;
     493        double  vz;
     494        double  maxvel;
     495
     496        /*retrive velocity values at nodes */
     497        inputs->GetParameterValue(&vx,VxEnum);
     498        inputs->GetParameterValue(&vy,VyEnum);
     499        if(dim==3)inputs->GetParameterValue(&vz,VzEnum);
     500
     501        /*now, compute maximum of velocity :*/
     502        if(dim==2){
     503                maxvel=sqrt(pow(vx,2)+pow(vy,2));
     504        }
     505        else{
     506                maxvel=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
     507        }
     508
     509        /*Assign output pointers:*/
     510        *pmaxvel=maxvel;
     511
     512}
     513/*}}}*/
     514/*FUNCTION Sing::MaxVx(double* pmaxvx, bool process_units);{{{1*/
     515void  Sing::MaxVx(double* pmaxvx, bool process_units){
     516
     517        int dim;
     518        double  maxvx;
     519
     520        /*retrieve dim parameter: */
     521        parameters->FindParam(&dim,DimEnum);
     522
     523        /*retrive velocity values at nodes */
     524        inputs->GetParameterValue(&maxvx,VxEnum);
     525
     526        /*Assign output pointers:*/
     527        *pmaxvx=maxvx;
     528
     529}
     530/*}}}*/
     531/*FUNCTION Sing::MaxVy(double* pmaxvy, bool process_units);{{{1*/
     532void  Sing::MaxVy(double* pmaxvy, bool process_units){
     533
     534        int dim;
     535        double  maxvy;
     536
     537        /*retrieve dim parameter: */
     538        parameters->FindParam(&dim,DimEnum);
     539
     540        /*retrive velocity values at nodes */
     541        inputs->GetParameterValue(&maxvy,VyEnum);
     542
     543        /*Assign output pointers:*/
     544        *pmaxvy=maxvy;
     545
     546}
     547/*}}}*/
     548/*FUNCTION Sing::MaxVz(double* pmaxvz, bool process_units);{{{1*/
     549void  Sing::MaxVz(double* pmaxvz, bool process_units){
     550
     551        int dim;
     552        double  maxvz;
     553
     554        /*retrieve dim parameter: */
     555        parameters->FindParam(&dim,DimEnum);
     556
     557        /*retrive velocity values at nodes */
     558        inputs->GetParameterValue(&maxvz,VzEnum);
     559
     560        /*Assign output pointers:*/
     561        *pmaxvz=maxvz;
     562
     563}
     564/*}}}*/
     565/*FUNCTION Sing::MinVel(double* pminvel, bool process_units);{{{1*/
     566void  Sing::MinVel(double* pminvel, bool process_units){
     567
     568        int dim;
     569        double  vx;
     570        double  vy;
     571        double  vz;
     572        double  minvel;
     573
     574        /*retrive velocity values at nodes */
     575        inputs->GetParameterValue(&vx,VxEnum);
     576        inputs->GetParameterValue(&vy,VyEnum);
     577        if(dim==3)inputs->GetParameterValue(&vz,VzEnum);
     578
     579        /*now, compute minimum of velocity :*/
     580        if(dim==2){
     581                minvel=sqrt(pow(vx,2)+pow(vy,2));
     582        }
     583        else{
     584                minvel=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
     585        }
     586
     587        /*Assign output pointers:*/
     588        *pminvel=minvel;
     589
     590}
     591/*}}}*/
     592/*FUNCTION Sing::MinVx(double* pminvx, bool process_units);{{{1*/
     593void  Sing::MinVx(double* pminvx, bool process_units){
     594
     595        int dim;
     596        double  minvx;
     597
     598        /*retrieve dim parameter: */
     599        parameters->FindParam(&dim,DimEnum);
     600
     601        /*retrive velocity values at nodes */
     602        inputs->GetParameterValue(&minvx,VxEnum);
     603
     604        /*Assign output pointers:*/
     605        *pminvx=minvx;
     606
     607}
     608/*}}}*/
     609/*FUNCTION Sing::MinVy(double* pminvy, bool process_units);{{{1*/
     610void  Sing::MinVy(double* pminvy, bool process_units){
     611
     612        int dim;
     613        double  minvy;
     614
     615        /*retrieve dim parameter: */
     616        parameters->FindParam(&dim,DimEnum);
     617
     618        /*retrive velocity values at nodes */
     619        inputs->GetParameterValue(&minvy,VyEnum);
     620
     621        /*Assign output pointers:*/
     622        *pminvy=minvy;
     623
     624}
     625/*}}}*/
     626/*FUNCTION Sing::MinVz(double* pminvz, bool process_units);{{{1*/
     627void  Sing::MinVz(double* pminvz, bool process_units){
     628
     629        int dim;
     630        double  minvz;
     631
     632        /*retrieve dim parameter: */
     633        parameters->FindParam(&dim,DimEnum);
     634
     635        /*retrive velocity values at nodes */
     636        inputs->GetParameterValue(&minvz,VzEnum);
     637
     638        /*Assign output pointers:*/
     639        *pminvz=minvz;
     640
     641}
     642/*}}}*/
     643/*FUNCTION Sing::Misfit {{{1*/
     644double Sing::Misfit(void){
     645        ISSMERROR(" not supported yet!");
     646}
     647/*}}}*/
     648/*FUNCTION Sing::PatchFill(int* pcount, Patch* patch){{{1*/
     649void  Sing::PatchFill(int* pcount, Patch* patch){
     650       
     651        ISSMERROR(" not supported yet!");
     652}
     653/*}}}*/
     654/*FUNCTION Sing::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){{{1*/
     655void  Sing::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){
     656
     657        ISSMERROR(" not supported yet!");
     658       
     659}
     660/*}}}*/
     661/*FUNCTION Sing::ProcessResultsUnits(void){{{1*/
     662void  Sing::ProcessResultsUnits(void){
     663        ISSMERROR(" not supported yet!");
     664}
     665/*}}}*/
     666/*FUNCTION Sing::SurfaceArea {{{1*/
     667double Sing::SurfaceArea( void){
     668        ISSMERROR(" not supported yet!");
     669}
     670/*}}}*/
     671/*FUNCTION Sing::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);{{{1*/
     672void Sing::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){
     673        ISSMERROR(" not supported yet!");
     674}
     675/*}}}*/
     676
     677/*Sing specific routines: */
    138678/*FUNCTION Sing::IsInput{{{1*/
    139679bool Sing::IsInput(int name){
    … …  
    143683        }
    144684        else return false;
    145 }
    146 /*}}}*/
    147 /*FUNCTION Sing::InputUpdateFromSolution {{{1*/
    148 void  Sing::InputUpdateFromSolution(double* solution){
    149         ISSMERROR(" not supported yet!");
    150 }
    151 /*}}}*/
    152 /*FUNCTION Sing::GetSolutionFromInputs(Vec solution);{{{1*/
    153 void  Sing::GetSolutionFromInputs(Vec solution){
    154         ISSMERROR(" not supported yet!");
    155 }
    156 /*}}}*/
    157 /*FUNCTION Sing::InputToResult(int enum_type,int step,double time){{{1*/
    158 void  Sing::InputToResult(int enum_type,int step,double time){
    159         ISSMERROR(" not supported yet!");
    160 }
    161 /*}}}*/
    162 /*FUNCTION Sing::ProcessResultsUnits(void){{{1*/
    163 void  Sing::ProcessResultsUnits(void){
    164         ISSMERROR(" not supported yet!");
    165 }
    166 /*}}}*/
    167                
    168 /*Sing functions*/
    169 /*FUNCTION Sing::ComputeBasalStress {{{1*/
    170 void  Sing::ComputeBasalStress(Vec p_g){
    171 
    172         ISSMERROR("Not implemented yet");
    173 
    174 }
    175 /*}}}*/
    176 /*FUNCTION Sing::ComputePressure {{{1*/
    177 void  Sing::ComputePressure(Vec p_g){
    178 
    179         int    dof;
    180         double pressure;
    181         double thickness;
    182         double rho_ice,g;
    183 
    184         /*Get dof list on which we will plug the pressure values: */
    185         GetDofList1(&dof);
    186 
    187         /*pressure is lithostatic: */
    188         rho_ice=matpar->GetRhoIce();
    189         g=matpar->GetG();
    190         inputs->GetParameterValue(&thickness,ThicknessEnum);
    191         pressure=rho_ice*g*thickness;
    192        
    193         /*plug local pressure values into global pressure vector: */
    194         VecSetValue(p_g,dof,pressure,INSERT_VALUES);
    195 
    196 }
    197 /*}}}*/
    198 /*FUNCTION Sing::ComputeStrainRate {{{1*/
    199 void  Sing::ComputeStrainRate(Vec p_g){
    200 
    201         ISSMERROR("Not implemented yet");
    202 
    203 }
    204 /*}}}*/
    205 /*FUNCTION Sing::CostFunction {{{1*/
    206 double Sing::CostFunction(){
    207         ISSMERROR(" not supported yet!");
    208 }
    209 /*}}}*/
    210 /*FUNCTION Sing::CreateKMatrix {{{1*/
    211 
    212 void  Sing::CreateKMatrix(Mat Kgg){
    213 
    214         int analysis_type;
    215 
    216         /*retrive parameters: */
    217         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    218 
    219         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    220         if (analysis_type==DiagnosticHutterAnalysisEnum){
    221                 CreateKMatrixDiagnosticHutter( Kgg);
    222         }
    223         else{
    224                 ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
    225         }
    226 
    227685}
    228686/*}}}*/
    … …  
    246704
    247705        MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
    248 
    249 }
    250 /*}}}*/
    251 /*FUNCTION Sing::CreatePVector {{{1*/
    252 void  Sing::CreatePVector(Vec pg){
    253 
    254         int analysis_type;
    255 
    256         /*retrive parameters: */
    257         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    258        
    259         /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
    260         if (analysis_type==DiagnosticHutterAnalysisEnum){
    261                         CreatePVectorDiagnosticHutter( pg);
    262         }
    263         else{
    264                 ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
    265         }
    266706
    267707}
    … …  
    317757}
    318758/*}}}*/
    319 /*FUNCTION Sing::Du {{{1*/
    320 void  Sing::Du(Vec){
    321         ISSMERROR(" not supported yet!");
    322 }
    323 /*}}}*/
    324 /*FUNCTION Sing::GetBedList {{{1*/
    325 void  Sing::GetBedList(double*){
    326         ISSMERROR(" not supported yet!");
    327 }
    328 /*}}}*/
    329759/*FUNCTION Sing::GetDofList {{{1*/
    330760void  Sing::GetDofList(int* doflist,int* pnumberofdofspernode){
    … …  
    352782}
    353783/*}}}*/
    354 /*FUNCTION Sing::GetMatPar {{{1*/
    355 void* Sing::GetMatPar(){
    356 
    357         return matpar;
    358 }
    359 /*}}}*/
    360 /*FUNCTION Sing::GetNodes {{{1*/
    361 void  Sing::GetNodes(void** vpnodes){
    362        
    363         Node** pnodes=NULL;
    364 
    365         /*recover nodes: */
    366         pnodes=(Node**)vpnodes;
    367 
    368         pnodes[0]=node;
    369 }
    370 /*}}}*/
    371 /*FUNCTION Sing::GetOnBed {{{1*/
    372 bool   Sing::GetOnBed(){
    373         ISSMERROR(" not supported yet!");
    374 }
    375 /*}}}*/
    376 /*FUNCTION Sing::GetShelf {{{1*/
    377 bool   Sing::GetShelf(){
    378         ISSMERROR(" not supported yet!");
    379 }
    380 /*}}}*/
    381 /*FUNCTION Sing::GetThicknessList {{{1*/
    382 void  Sing::GetThicknessList(double* thickness_list){
    383         ISSMERROR(" not supported yet!");
    384 }
    385 /*}}}*/
    386 /*FUNCTION Sing::Gradj {{{1*/
    387 void  Sing::Gradj(Vec gradient,int control_type){
    388         ISSMERROR(" not supported yet!");
    389 }
    390 /*}}}*/
    391 /*FUNCTION Sing::GradB {{{1*/
    392 void  Sing::GradjB(Vec gradient){
    393         ISSMERROR(" not supported yet!");
    394 }
    395 /*}}}*/
    396 /*FUNCTION Sing::GradjDrag {{{1*/
    397 void  Sing::GradjDrag(Vec gradient){
    398         ISSMERROR(" not supported yet!");
    399 }
    400 /*}}}*/
    401 /*FUNCTION Sing::MassFlux {{{1*/
    402 double Sing::MassFlux( double* segment){
    403         ISSMERROR(" not supported yet!");
    404 }
    405 /*}}}*/
    406 /*FUNCTION Sing::Misfit {{{1*/
    407 double Sing::Misfit(void){
    408         ISSMERROR(" not supported yet!");
    409 }
    410 /*}}}*/
    411 /*FUNCTION Sing::SurfaceArea {{{1*/
    412 double Sing::SurfaceArea( void){
    413         ISSMERROR(" not supported yet!");
    414 }
    415 /*}}}*/
    416784/*FUNCTION Sing::SetClone {{{1*/
    417785void  Sing::SetClone(int* minranks){
    … …  
    420788}
    421789/*}}}1*/
    422 /*FUNCTION Sing::InputUpdateFromVector(double* vector, int name, int type);{{{1*/
    423 void  Sing::InputUpdateFromVector(double* vector, int name, int type){
    424 
    425         /*Check that name is an element input*/
    426         if (!IsInput(name)) return;
    427 
    428         switch(type){
    429 
    430                 case VertexEnum:
    431 
    432                         /*New SingVertexInpu*/
    433                         double value;
    434 
    435                         /*Get values on the 6 vertices*/
    436                         value=vector[node->GetVertexDof()];
    437 
    438                         /*update input*/
    439                         this->inputs->AddInput(new SingVertexInput(name,value));
    440                         return;
    441 
    442                 default:
    443                         ISSMERROR("type %i (%s) not implemented yet",type,EnumAsString(type));
    444         }
    445 }
    446 /*}}}*/
    447 /*FUNCTION Sing::InputUpdateFromVector(int* vector, int name, int type);{{{1*/
    448 void  Sing::InputUpdateFromVector(int* vector, int name, int type){
    449         ISSMERROR(" not supported yet!");
    450 }
    451 /*}}}*/
    452 /*FUNCTION Sing::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/
    453 void  Sing::InputUpdateFromVector(bool* vector, int name, int type){
    454         ISSMERROR(" not supported yet!");
    455 }
    456 /*}}}*/
    457 /*FUNCTION Sing::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){{{1*/
    458 void  Sing::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){
    459 
    460         ISSMERROR(" not supported yet!");
    461        
    462 }
    463 /*}}}*/
    464 /*FUNCTION Sing::PatchFill(int* pcount, Patch* patch){{{1*/
    465 void  Sing::PatchFill(int* pcount, Patch* patch){
    466        
    467         ISSMERROR(" not supported yet!");
    468 }
    469 /*}}}*/
    470 /*FUNCTION Sing::MinVel(double* pminvel, bool process_units);{{{1*/
    471 void  Sing::MinVel(double* pminvel, bool process_units){
    472 
    473         int dim;
    474         double  vx;
    475         double  vy;
    476         double  vz;
    477         double  minvel;
    478 
    479         /*retrive velocity values at nodes */
    480         inputs->GetParameterValue(&vx,VxEnum);
    481         inputs->GetParameterValue(&vy,VyEnum);
    482         if(dim==3)inputs->GetParameterValue(&vz,VzEnum);
    483 
    484         /*now, compute minimum of velocity :*/
    485         if(dim==2){
    486                 minvel=sqrt(pow(vx,2)+pow(vy,2));
    487         }
    488         else{
    489                 minvel=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
    490         }
    491 
    492         /*Assign output pointers:*/
    493         *pminvel=minvel;
    494 
    495 }
    496 /*}}}*/
    497 /*FUNCTION Sing::MaxVel(double* pmaxvel, bool process_units);{{{1*/
    498 void  Sing::MaxVel(double* pmaxvel, bool process_units){
    499 
    500         int dim;
    501         double  vx;
    502         double  vy;
    503         double  vz;
    504         double  maxvel;
    505 
    506         /*retrive velocity values at nodes */
    507         inputs->GetParameterValue(&vx,VxEnum);
    508         inputs->GetParameterValue(&vy,VyEnum);
    509         if(dim==3)inputs->GetParameterValue(&vz,VzEnum);
    510 
    511         /*now, compute maximum of velocity :*/
    512         if(dim==2){
    513                 maxvel=sqrt(pow(vx,2)+pow(vy,2));
    514         }
    515         else{
    516                 maxvel=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
    517         }
    518 
    519         /*Assign output pointers:*/
    520         *pmaxvel=maxvel;
    521 
    522 }
    523 /*}}}*/
    524 /*FUNCTION Sing::MinVx(double* pminvx, bool process_units);{{{1*/
    525 void  Sing::MinVx(double* pminvx, bool process_units){
    526 
    527         int dim;
    528         double  minvx;
    529 
    530         /*retrieve dim parameter: */
    531         parameters->FindParam(&dim,DimEnum);
    532 
    533         /*retrive velocity values at nodes */
    534         inputs->GetParameterValue(&minvx,VxEnum);
    535 
    536         /*Assign output pointers:*/
    537         *pminvx=minvx;
    538 
    539 }
    540 /*}}}*/
    541 /*FUNCTION Sing::MaxVx(double* pmaxvx, bool process_units);{{{1*/
    542 void  Sing::MaxVx(double* pmaxvx, bool process_units){
    543 
    544         int dim;
    545         double  maxvx;
    546 
    547         /*retrieve dim parameter: */
    548         parameters->FindParam(&dim,DimEnum);
    549 
    550         /*retrive velocity values at nodes */
    551         inputs->GetParameterValue(&maxvx,VxEnum);
    552 
    553         /*Assign output pointers:*/
    554         *pmaxvx=maxvx;
    555 
    556 }
    557 /*}}}*/
    558 /*FUNCTION Sing::MaxAbsVx(double* pmaxabsvx, bool process_units);{{{1*/
    559 void  Sing::MaxAbsVx(double* pmaxabsvx, bool process_units){
    560 
    561         int dim;
    562         double  maxabsvx;
    563 
    564         /*retrieve dim parameter: */
    565         parameters->FindParam(&dim,DimEnum);
    566 
    567         /*retrive velocity values at nodes */
    568         inputs->GetParameterValue(&maxabsvx,VxEnum);
    569         maxabsvx=fabs(maxabsvx);
    570 
    571         /*Assign output pointers:*/
    572         *pmaxabsvx=maxabsvx;
    573 }
    574 /*}}}*/
    575 /*FUNCTION Sing::MinVy(double* pminvy, bool process_units);{{{1*/
    576 void  Sing::MinVy(double* pminvy, bool process_units){
    577 
    578         int dim;
    579         double  minvy;
    580 
    581         /*retrieve dim parameter: */
    582         parameters->FindParam(&dim,DimEnum);
    583 
    584         /*retrive velocity values at nodes */
    585         inputs->GetParameterValue(&minvy,VyEnum);
    586 
    587         /*Assign output pointers:*/
    588         *pminvy=minvy;
    589 
    590 }
    591 /*}}}*/
    592 /*FUNCTION Sing::MaxVy(double* pmaxvy, bool process_units);{{{1*/
    593 void  Sing::MaxVy(double* pmaxvy, bool process_units){
    594 
    595         int dim;
    596         double  maxvy;
    597 
    598         /*retrieve dim parameter: */
    599         parameters->FindParam(&dim,DimEnum);
    600 
    601         /*retrive velocity values at nodes */
    602         inputs->GetParameterValue(&maxvy,VyEnum);
    603 
    604         /*Assign output pointers:*/
    605         *pmaxvy=maxvy;
    606 
    607 }
    608 /*}}}*/
    609 /*FUNCTION Sing::MaxAbsVy(double* pmaxabsvy, bool process_units);{{{1*/
    610 void  Sing::MaxAbsVy(double* pmaxabsvy, bool process_units){
    611 
    612         int dim;
    613         double  maxabsvy;
    614 
    615         /*retrieve dim parameter: */
    616         parameters->FindParam(&dim,DimEnum);
    617 
    618         /*retrive velocity values at nodes */
    619         inputs->GetParameterValue(&maxabsvy,VyEnum);
    620         maxabsvy=fabs(maxabsvy);
    621 
    622         /*Assign output pointers:*/
    623         *pmaxabsvy=maxabsvy;
    624 }
    625 /*}}}*/
    626 /*FUNCTION Sing::MinVz(double* pminvz, bool process_units);{{{1*/
    627 void  Sing::MinVz(double* pminvz, bool process_units){
    628 
    629         int dim;
    630         double  minvz;
    631 
    632         /*retrieve dim parameter: */
    633         parameters->FindParam(&dim,DimEnum);
    634 
    635         /*retrive velocity values at nodes */
    636         inputs->GetParameterValue(&minvz,VzEnum);
    637 
    638         /*Assign output pointers:*/
    639         *pminvz=minvz;
    640 
    641 }
    642 /*}}}*/
    643 /*FUNCTION Sing::MaxVz(double* pmaxvz, bool process_units);{{{1*/
    644 void  Sing::MaxVz(double* pmaxvz, bool process_units){
    645 
    646         int dim;
    647         double  maxvz;
    648 
    649         /*retrieve dim parameter: */
    650         parameters->FindParam(&dim,DimEnum);
    651 
    652         /*retrive velocity values at nodes */
    653         inputs->GetParameterValue(&maxvz,VzEnum);
    654 
    655         /*Assign output pointers:*/
    656         *pmaxvz=maxvz;
    657 
    658 }
    659 /*}}}*/
    660 /*FUNCTION Sing::MaxAbsVz(double* pmaxabsvz, bool process_units);{{{1*/
    661 void  Sing::MaxAbsVz(double* pmaxabsvz, bool process_units){
    662 
    663         int dim;
    664         double  maxabsvz;
    665 
    666         /*retrieve dim parameter: */
    667         parameters->FindParam(&dim,DimEnum);
    668 
    669         /*retrive velocity values at nodes */
    670         inputs->GetParameterValue(&maxabsvz,VzEnum);
    671         maxabsvz=fabs(maxabsvz);
    672 
    673         /*Assign output pointers:*/
    674         *pmaxabsvz=maxabsvz;
    675 }
    676 /*}}}*/
    677 /*FUNCTION Sing::InputDuplicate(int original_enum,int new_enum){{{1*/
    678 void  Sing::InputDuplicate(int original_enum,int new_enum){
    679 
    680         Input* original=NULL;
    681         Input* copy=NULL;
    682 
    683         /*Make a copy of the original input: */
    684         original=(Input*)this->inputs->GetInput(original_enum);
    685         copy=(Input*)original->copy();
    686 
    687         /*Change copy enum to reinitialized_enum: */
    688         copy->ChangeEnum(new_enum);
    689 
    690         /*Add copy into inputs, it will wipe off the one already there: */
    691         inputs->AddObject((Input*)copy);
    692 }
    693 /*}}}*/
    694 /*FUNCTION Sing::InputScale(int enum_type,double scale_factor){{{1*/
    695 void  Sing::InputScale(int enum_type,double scale_factor){
    696 
    697         Input* input=NULL;
    698 
    699         /*Make a copy of the original input: */
    700         input=(Input*)this->inputs->GetInput(enum_type);
    701 
    702         /*Scale: */
    703         input->Scale(scale_factor);
    704 }
    705 /*}}}*/
    706 /*FUNCTION Sing::InputAXPY(int YEnum, double scalar, int XEnum);{{{1*/
    707 void  Sing::InputAXPY(int YEnum, double scalar, int XEnum){
    708 
    709         Input* xinput=NULL;
    710         Input* yinput=NULL;
    711 
    712         /*Find x and y inputs: */
    713         xinput=(Input*)this->inputs->GetInput(XEnum);
    714         yinput=(Input*)this->inputs->GetInput(YEnum);
    715 
    716         /*some checks: */
    717         if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!");
    718         if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!");
    719 
    720         /*Scale: */
    721         yinput->AXPY(xinput,scalar);
    722 }
    723 /*}}}*/
    724 /*FUNCTION Sing::InputControlConstrain(int control_type, double cm_min, double cm_max){{{1*/
    725 void  Sing::InputControlConstrain(int control_type, double cm_min, double cm_max){
    726 
    727         Input* input=NULL;
    728 
    729         /*Find input: */
    730         input=(Input*)this->inputs->GetInput(control_type);
    731        
    732         /*Do nothing if we  don't find it: */
    733         if(!input)return;
    734 
    735         /*Constrain input using cm_min and cm_max: */
    736         input->Constrain(cm_min,cm_max);
    737 
    738 }
    739 /*}}}*/
    740 /*FUNCTION Sing::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/
    741 void  Sing::GetVectorFromInputs(Vec vector,int NameEnum){
    742 
    743         int i;
    744         const int numvertices=1;
    745         int doflist1[numvertices];
    746 
    747         /*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
    748         for(i=0;i<this->inputs->Size();i++){
    749                 Input* input=(Input*)this->inputs->GetObjectByOffset(i);
    750                 if(input->EnumType()==NameEnum){
    751                         /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
    752                         this->GetDofList1(&doflist1[0]);
    753                         input->GetVectorFromInputs(vector,&doflist1[0]);
    754                         break;
    755                 }
    756         }
    757 }
    758 /*}}}*/
    759 /*FUNCTION Sing::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/
    760 void  Sing::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){
    761 
    762         int i;
    763         Input** new_inputs=NULL;
    764         Input** old_inputs=NULL;
    765         int     converged=1;
    766 
    767         new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
    768         old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
    769        
    770         for(i=0;i<num_enums/2;i++){
    771                 new_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+0]);
    772                 old_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+1]);
    773                 if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(enums[2*i+0]));
    774                 if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(enums[2*i+0]));
    775         }
    776 
    777         /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
    778         for(i=0;i<num_criterionenums;i++){
    779                 IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,criterionenums[i]);
    780                 if(eps[i]>criterionvalues[i]) converged=0;
    781         }
    782 
    783         /*Assign output pointers:*/
    784         *pconverged=converged;
    785 
    786 }
    787 /*}}}*/
Note: See TracChangeset for help on using the changeset viewer.