Changeset 1805


Ignore:
Timestamp:
08/24/09 14:26:14 (16 years ago)
Author:
Eric.Larour
Message:

New rift non-linear convergence schema implemented-> increase fraction until convergence is reached. src/c/objects/Riftfront.cpp src/c/objects/Riftfront.h src/c/PenaltyConstraintsx/RiftConstraints.cpp
New parallel thermal static and thermalstatic control solutions. src/c/parallel/thermalstatic.cpp src/c/parallel/thermalstatic_core.cpp src/c/parallel/objectivefunctionC.cpp src/c/parallel/GradJCompute.cpp src/c/parallel/ControlInitialization.cpp src/c/parallel/diagnostic_core.cpp
New OutputRifts module to output penalty fractions. src/c/OutputRiftsx/OutputRiftsx.cpp src/c/OutputRiftsx/OutputRiftsx.h src/c/Makefile.am src/mex/OutputRifts src/mex/OutputRifts/OutputRifts.cpp src/mex/OutputRifts/OutputRifts.h src/mex/TriMeshProcessRifts/TriMeshProcessRifts.cpp src/mex/Makefile.am

Location:
issm/trunk/src
Files:
9 added
29 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/DataSet/DataSet.cpp

    r1784 r1805  
    13831383
    13841384}
     1385               
     1386void  DataSet::OutputRifts(Vec riftproperties){
     1387
     1388        vector<Object*>::iterator object;
     1389        Riftfront* riftfront=NULL;
     1390
     1391        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1392
     1393                if((*object)->Enum()==RiftfrontEnum()){
     1394
     1395                        riftfront=(Riftfront*)(*object);
     1396                        riftfront->OutputProperties(riftproperties);
     1397                }
     1398        }
     1399
     1400
     1401
     1402}

  • issm/trunk/src/c/DataSet/DataSet.h

    r1784 r1805  
    8282                void  FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse);
    8383                void  UpdateNodePositions(double* thickness,double* bed);
     84                void  OutputRifts(Vec riftproperties);
    8485
    8586};

  • issm/trunk/src/c/Makefile.am

    r1787 r1805  
    292292                                        ./FieldDepthAveragex/FieldDepthAveragex.h\
    293293                                        ./FieldExtrudex/FieldExtrudex.cpp\
    294                                         ./FieldExtrudex/FieldExtrudex.h
     294                                        ./FieldExtrudex/FieldExtrudex.h\
     295                                        ./OutputRiftsx/OutputRiftsx.h\
     296                                        ./OutputRiftsx/OutputRiftsx.cpp
    295297
    296298
    … …  
    577579                                        ./FieldExtrudex/FieldExtrudex.cpp\
    578580                                        ./FieldExtrudex/FieldExtrudex.h\
     581                                        ./OutputRiftsx/OutputRiftsx.h\
     582                                        ./OutputRiftsx/OutputRiftsx.cpp\
    579583                                        ./parallel/diagnostic_core.cpp\
    580584                                        ./parallel/diagnostic_core_linear.cpp\
    … …  
    596600                                        ./parallel/transient_core_2d.cpp\
    597601                                        ./parallel/transient_core_3d.cpp\
     602                                        ./parallel/thermalstatic_core.cpp\
    598603                                        ./parallel/OutputResults.cpp\
    599604                                        ./parallel/OutputControl.cpp
    … …  
    604609bin_PROGRAMS =
    605610else
    606 bin_PROGRAMS = diagnostic.exe  control.exe thermal.exe prognostic.exe transient.exe
     611bin_PROGRAMS = diagnostic.exe  control.exe thermal.exe prognostic.exe transient.exe thermalstatic.exe
    607612endif
    608613
    … …  
    612617diagnostic_exe_CXXFLAGS= -fPIC -D_PARALLEL_
    613618
     619thermalstatic_exe_SOURCES = parallel/thermalstatic.cpp
     620thermalstatic_exe_CXXFLAGS= -fPIC -D_PARALLEL_
     621
    614622control_exe_SOURCES = parallel/control.cpp
    615623control_exe_CXXFLAGS= -fPIC -D_PARALLEL_

  • issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp

    r1648 r1805  
    1818
    1919void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,Model* model,ConstDataHandle model_handle){
    20 
    2120
    2221        /*create parameters common to all solutions: */
    … …  
    125124                throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s"," analysis_type: ",model->analysis_type," sub_analysis_type: ",model->sub_analysis_type," not supported yet!"));
    126125        }
     126
    127127                       
    128128}

  • issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp

    r1666 r1805  
    203203        param->SetInteger(numberofdofspernode);
    204204        parameters->AddObject(param);
     205
     206        /*numrifts: */
     207        ModelFetchData((void**)&model->riftinfo,&model->numrifts,NULL,model_handle,"riftinfo","Matrix","Mat");
     208        count++;
     209        param= new Param(count,"numrifts",INTEGER);
     210        param->SetInteger(model->numrifts);
     211        parameters->AddObject(param);
     212        xfree((void**)&model->riftinfo);
    205213       
    206214        /*All our datasets are already ordered by ids. Set presort flag so that later on, when sorting is requested on these

  • issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp

    r1784 r1805  
    5555        double riftfront_friction;
    5656        double riftfront_fraction;
     57        double riftfront_fractionincrement;
    5758        bool   riftfront_shelf;
    5859        double riftfront_penalty_offset;
    … …  
    6869        double friction;
    6970        double fraction;
     71        double fractionincrement;
    7072       
    7173        /*Create loads: */
    … …  
    186188                        friction=*(model->riftinfo+RIFTINFOSIZE*i+8);
    187189                        fraction=*(model->riftinfo+RIFTINFOSIZE*i+9);
     190                        fractionincrement=*(model->riftinfo+RIFTINFOSIZE*i+10);
    188191       
    189192                        strcpy(riftfront_type,"2d");
    … …  
    209212                        riftfront_friction=friction;
    210213                        riftfront_fraction=fraction;
     214                        riftfront_fractionincrement=fractionincrement;
    211215                        riftfront_shelf=(bool)model->gridoniceshelf[grid1-1];
    212216
    … …  
    218222                        riftfront_prestable=0;
    219223                       
    220                         riftfront=new Riftfront(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_counter,riftfront_prestable,riftfront_shelf);
     224                        riftfront=new Riftfront(riftfront_type,riftfront_id, riftfront_node_ids, riftfront_mparid, riftfront_h,riftfront_b,riftfront_s,riftfront_normal,riftfront_length,riftfront_fill,riftfront_friction, riftfront_fraction, riftfront_fractionincrement, riftfront_penalty_offset, riftfront_penalty_lock, riftfront_active,riftfront_counter,riftfront_prestable,riftfront_shelf);
    221225
    222226                        loads->AddObject(riftfront);

  • issm/trunk/src/c/ModelProcessorx/Model.h

    r1784 r1805  
    1111#include "../toolkits/toolkits.h"
    1212
    13 #define RIFTINFOSIZE 10
     13#define RIFTINFOSIZE 11
    1414
    1515struct Model {

  • issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp

    r1784 r1805  
    1616        int potential;
    1717
    18         /*First, has the material non-linearity stabilized? : */
    19         if(!IsMaterialStable(loads,inputs,analysis_type)){
    20                 /*Do nothing, no constraints activated!: */
    21                 printf("converging material\n");
    22                 num_unstable_constraints=0;
    23                 converged=0;
    24         }
    25         else{
    26                 /*Ok, start constraining: */
    27                 printf("converged material: constraining\n");
    28 
    29                 #ifdef _ZIGZAGCOUNTER_
    30                 if(!IsPreStable(loads)){
    31 
    32                         PreConstrain(&num_unstable_constraints,loads,inputs,analysis_type);
    33 
    34                         if(!num_unstable_constraints){
    35                                 /*Ok, preconstraining is done. Figure out size of unstable constraints pool. Then set
    36                                  * prestable flag to 1 for all rift grids: */
    37 
    38                                 potential=PotentialUnstableConstraints(loads,inputs,analysis_type);
    39                                 printf("      set of constraints is prestabilised; unstable constraints potentialnumber: %i\n",potential);
    40 
    41                                 SetPreStable(loads);
    42                                 num_unstable_constraints=-1; //ensure convergence does not happen!
    43                         }
    44 
    45                 }
    46                 else{
    47                         /*Ok, set of constraints is stable. Start increasing number of inactive constraints, one by one: */
    48                
    49                         /*Figure out max penetration: */
    50                         MaxPenetrationInInputs(loads,inputs,analysis_type);
    51                        
    52                         /*Constrain rifts: */
    53                         Constrain(&num_unstable_constraints,loads,inputs,analysis_type);
    54                 }
    55                 #else
    56                         Constrain(&num_unstable_constraints,loads,inputs,analysis_type);
    57                 #endif
    58 
    59                 if(num_unstable_constraints==0)converged=1;
    60         }
     18        Constrain(&num_unstable_constraints,loads,inputs,analysis_type);
     19        if(num_unstable_constraints==0)converged=1;
    6120
    6221        /*Assign output pointers: */

  • issm/trunk/src/c/issm.h

    r1665 r1805  
    5757#include "./NodeConnectivityx/NodeConnectivityx.h"
    5858#include "./ElementConnectivityx/ElementConnectivityx.h"
     59#include "./OutputRiftsx/OutputRiftsx.h"
    5960
    6061

  • issm/trunk/src/c/objects/Riftfront.cpp

    r1784 r1805  
    2525}
    2626
    27 Riftfront::Riftfront(char riftfront_type[RIFTFRONTSTRING],int riftfront_id, int riftfront_node_ids[MAX_RIFTFRONT_GRIDS], int riftfront_mparid, double riftfront_h[MAX_RIFTFRONT_GRIDS],double riftfront_b[MAX_RIFTFRONT_GRIDS],double riftfront_s[MAX_RIFTFRONT_GRIDS],double riftfront_normal[2],double riftfront_length,int riftfront_fill,double riftfront_friction, double riftfront_fraction,double riftfront_penalty_offset, int riftfront_penalty_lock, bool riftfront_active,int riftfront_counter,bool riftfront_prestable,bool riftfront_shelf){
     27Riftfront::Riftfront(char riftfront_type[RIFTFRONTSTRING],int riftfront_id, int riftfront_node_ids[MAX_RIFTFRONT_GRIDS], int riftfront_mparid, double riftfront_h[MAX_RIFTFRONT_GRIDS],double riftfront_b[MAX_RIFTFRONT_GRIDS],double riftfront_s[MAX_RIFTFRONT_GRIDS],double riftfront_normal[2],double riftfront_length,int riftfront_fill,double riftfront_friction, double riftfront_fraction,double riftfront_fractionincrement, double riftfront_penalty_offset, int riftfront_penalty_lock, bool riftfront_active,int riftfront_counter,bool riftfront_prestable,bool riftfront_shelf){
    2828
    2929        int i;
    … …  
    5454        friction=riftfront_friction;
    5555        fraction=riftfront_fraction;
     56        fractionincrement=riftfront_fractionincrement;
    5657        penalty_offset=riftfront_penalty_offset;
    5758        penalty_lock=riftfront_penalty_lock;
    … …  
    8586
    8687        printf("normal [%g,%g], length %g\n",normal[0],normal[1],normal[2]);
    87         printf("fill: %i friction %g fraction %g\n",fill,friction,fraction);
     88        printf("fill: %i friction %g fraction %g fractionincrement %g \n",fill,friction,fraction,fractionincrement);
    8889        printf("penalty_offset %g\n",penalty_offset);
    8990        printf("penalty_lock %i\n",penalty_lock);
    … …  
    114115
    115116        printf("normal [%g,%g], length %g\n",normal[0],normal[1],normal[2]);
    116         printf("fill: %i friction %g fraction %g\n",fill,friction,fraction);
     117        printf("fill: %i friction %g fraction %g fractionincrement %g \n",fill,friction,fraction,fractionincrement);
    117118        printf("penalty_offset %g\n",penalty_offset);
    118119        printf("penalty_lock %i\n",penalty_lock);
    … …  
    156157        memcpy(marshalled_dataset,&friction,sizeof(friction));marshalled_dataset+=sizeof(friction);
    157158        memcpy(marshalled_dataset,&fraction,sizeof(fraction));marshalled_dataset+=sizeof(fraction);
     159        memcpy(marshalled_dataset,&fractionincrement,sizeof(fractionincrement));marshalled_dataset+=sizeof(fractionincrement);
    158160        memcpy(marshalled_dataset,&penalty_offset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
    159161        memcpy(marshalled_dataset,&penalty_lock,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock);
    … …  
    184186                sizeof(friction)+
    185187                sizeof(fraction)+
     188                sizeof(fractionincrement)+
    186189                sizeof(penalty_offset)+
    187190                sizeof(penalty_lock)+
    … …  
    225228        memcpy(&friction,marshalled_dataset,sizeof(friction));marshalled_dataset+=sizeof(friction);
    226229        memcpy(&fraction,marshalled_dataset,sizeof(fraction));marshalled_dataset+=sizeof(fraction);
     230        memcpy(&fractionincrement,marshalled_dataset,sizeof(fractionincrement));marshalled_dataset+=sizeof(fractionincrement);
    227231        memcpy(&penalty_offset,marshalled_dataset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
    228232        memcpy(&penalty_lock,marshalled_dataset,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock);
    … …  
    621625        penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1];
    622626
    623         #ifdef _ZIGZAGCOUNTER_
    624 
    625         found=inputs->Recover("max_penetration",&max_penetration);
    626         if(!found)throw ErrorException(__FUNCT__," could not find max_penetration in inputs!");
    627 
    628         /*Activate or deactivate penalties: */
    629         if(penetration<0){
    630                 //printf("riftfront %i penetrating.\n",this->GetId());
    631                 /*There is penetration, we need to active the penalty so this penetration will be NULL: */
    632                 activate=1;
    633         }
    634         else{
    635                 /*Ok, there is no penetration for these two grids of the rift. We want to deactivate  the penalty. If we do
    636                  * it all at once, then we will zigzag forever. Only deactivate once at a time. Deactivate the one that most
    637                  * wants to, ie the one with the max penetration: */
    638                 if (penetration==max_penetration){
    639                         //printf("riftfront %i at max penetration.\n",this->GetId());
    640                         activate=0;
    641                 }
    642                 else{
    643                         /*Ok, we are dealing with the  rest of the penalties that want to be freed. But may be in this lot there
    644                          * are penalties that were already free? Keep them as such. For the rest, do not release them yet: */
    645                         if (this->active==0){
    646                                 activate=0;
    647                         }
    648                         else{
    649                                 activate=1;
    650                         }
    651                 }
    652         }
    653 
    654         /*Here, we try to avoid zigzaging. When a penalty activates and deactivates for more than penalty_lock times,
    655          * we lock it: */
    656         if(this->counter>this->penalty_lock){
    657                 /*This penalty has zig zagged too long, fix it to smooth results: */
    658                 activate=1; this->active=activate;
    659                 if (this->counter==(this->penalty_lock+1)){
    660                         printf("riftfront %i locking at %g m/yr\n",this->GetId(),penetration*365.0*24.0*3600.0);
    661                         this->counter++;
    662                 }
    663         }
    664         #else
     627        /*activation: */
    665628        if(penetration<0)activate=1;
    666629        else  activate=0;
    667         #endif
     630
     631        /*Here, we try to avoid zigzaging. When a penalty activates and deactivates for more than penalty_lock times,
     632         * we increase the fraction of melange:*/
     633        if(this->counter>this->penalty_lock){
     634                /*reset counter: */
     635                this->counter=0;
     636                /*increase melange fraction: */
     637                this->fraction+=this->fractionincrement;
     638                if (this->fraction>1)this->fraction=(double)1.0;
     639                //printf("riftfront %i fraction: %g\n",this->GetId(),this->fraction);
     640        }
    668641
    669642        //Figure out stability of this penalty
    … …  
    814787        return this->material_converged;
    815788}
     789               
     790#undef __FUNCT__
     791#define __FUNCT__ "Riftfront::OutputProperties"
     792void  Riftfront::OutputProperties(Vec riftproperties){
     793
     794        int row_id=0;
     795        double value;
     796
     797        /*recover id of penalty: */
     798        row_id=this->GetId()-1; //c indexing, ids were matlab indexed
     799        value=(double)this->fraction;
     800
     801        /*Plug id and fraction  into riftproperties matrix: */
     802        VecSetValues(riftproperties,1,&row_id,&value,INSERT_VALUES);
     803}

  • issm/trunk/src/c/objects/Riftfront.h

    r1784 r1805  
    4242                double      friction;
    4343                double      fraction;
     44                double      fractionincrement;
    4445                bool        shelf;
    4546
    … …  
    5758
    5859                Riftfront();
    59                 Riftfront(char type[RIFTFRONTSTRING],int id, int node_ids[MAX_RIFTFRONT_GRIDS], int mparid, double h[MAX_RIFTFRONT_GRIDS],double b[MAX_RIFTFRONT_GRIDS],double s[MAX_RIFTFRONT_GRIDS],double normal[2],double length,int fill,double friction, double fraction, double penalty_offset, int penalty_lock,bool active,int counter,bool prestable,bool shelf);
     60                Riftfront(char type[RIFTFRONTSTRING],int id, int node_ids[MAX_RIFTFRONT_GRIDS], int mparid, double h[MAX_RIFTFRONT_GRIDS],double b[MAX_RIFTFRONT_GRIDS],double s[MAX_RIFTFRONT_GRIDS],double normal[2],double length,int fill,double friction, double fraction, double fractionincrement, double penalty_offset, int penalty_lock,bool active,int counter,bool prestable,bool shelf);
    6061                ~Riftfront();
    6162
    … …  
    8485                int   PotentialUnstableConstraint(int* punstable, void* inputs, int analysis_type);
    8586                int   IsMaterialStable(void* inputs, int analysis_type);
     87                void  OutputProperties(Vec riftproperties);
    8688                Object* copy();
    8789};

  • issm/trunk/src/c/parallel/ControlInitialization.cpp

    r1707 r1805  
    9797        //horizontal velocity
    9898        if(debug)_printf_("%s\n"," computing horizontal velocities...");
    99         diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
     99        diagnostic_core_nonlinear(&ug,NULL,NULL,NULL, fem_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
    100100        if(debug)_printf_("%s\n"," extruding horizontal velocities...");
    101101        VecDuplicatePatch(&ug_horiz,ug); FieldExtrudex( ug_horiz,fem_dh->elements,fem_dh->nodes, fem_dh->loads,fem_dh-> materials,"velocity",1);
    … …  
    138138        if(debug)_printf_("%s\n"," computing stokes velocities and pressure ...");
    139139        VecFree(&ug);
    140         diagnostic_core_nonlinear(&ug,NULL,NULL,fem_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
     140        diagnostic_core_nonlinear(&ug,NULL,NULL,NULL, fem_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
    141141        VecFree(&pg);
    142142

  • issm/trunk/src/c/parallel/GradJCompute.cpp

    r1210 r1805  
    5959
    6060        _printf_("%s\n","      recover solution for this stiffness and right hand side:");
    61         diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,femmodel,inputs,DiagnosticAnalysisEnum(),sub_analysis_type);
     61        diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,NULL, femmodel,inputs,DiagnosticAnalysisEnum(),sub_analysis_type);
    6262        VecToMPISerial(&u_g_double,u_g); VecFree(&u_g);
    6363        inputs->Add("velocity",u_g_double,numberofdofspernode,numberofnodes);

  • issm/trunk/src/c/parallel/ProcessResults.cpp

    r1209 r1805  
    7070        double* param_g=NULL;
    7171        double* parameter=NULL;
     72        Vec     riftproperties=NULL;
     73        double* riftproperties_serial=NULL;
     74        int     numrifts=0;
    7275
    7376        Vec     s_g=NULL;
    … …  
    106109                fem_dhu=fems+3;
    107110                fem_sl=fems+4;
     111       
     112                /*some flags needed: */
     113                fem_dh->parameters->FindParam((void*)&dim,"dim");
     114                fem_dhu->parameters->FindParam((void*)&ishutter,"ishutter");
     115                fem_ds->parameters->FindParam((void*)&isstokes,"isstokes");
     116                fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn");
     117        }
     118        if(analysis_type==ThermalstaticAnalysisEnum()){
     119                fem_dh=fems+0;
     120                fem_dv=fems+1;
     121                fem_ds=fems+2;
     122                fem_dhu=fems+3;
     123                fem_sl=fems+4;
     124                fem_t=fems+5;
    108125       
    109126                /*some flags needed: */
    … …  
    411428                        xfree((void**)&partition);
    412429                }
     430                else if(strcmp(result->GetFieldName(),"riftproperties")==0){
     431                        result->GetField(&riftproperties);
     432                        fem_dh->parameters->FindParam((void*)&numrifts,"numrifts");
     433                        VecToMPISerial(&riftproperties_serial,riftproperties);
     434                       
     435                        /*Ok, add parameter to newresults: */
     436                        newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"riftproperties",riftproperties_serial,numrifts);
     437                        newresults->AddObject(newresult);
     438
     439                }
    413440                else{
    414441                        /*Just copy the result into the new results dataset: */

  • issm/trunk/src/c/parallel/diagnostic_core.cpp

    r1723 r1805  
    3535        Vec slopex=NULL;
    3636        Vec slopey=NULL;
     37        Vec riftproperties=NULL;
    3738
    3839        /*flags: */
    … …  
    4647        int numberofdofspernode_ds;
    4748        int numberofnodes;
     49        int numrifts=0;
    4850
    4951        double stokesreconditioning;
    … …  
    7072        fem_ds->parameters->FindParam((void*)&isstokes,"isstokes");
    7173        fem_ds->parameters->FindParam((void*)&stokesreconditioning,"stokesreconditioning");
     74        fem_dh->parameters->FindParam((void*)&numrifts,"numrifts");
    7275
    7376        //specific parameters for specific models
    … …  
    112115               
    113116                if(debug)_printf_("%s\n"," computing horizontal velocities...");
    114                 diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
     117                diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh->loads,fem_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
    115118
    116119                if(debug)_printf_("%s\n"," computing pressure according to MacAyeal...");
    … …  
    169172                        if(debug)_printf_("%s\n"," computing stokes velocities and pressure ...");
    170173                        VecFree(&ug);
    171                         diagnostic_core_nonlinear(&ug,NULL,NULL,fem_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
     174                        diagnostic_core_nonlinear(&ug,NULL,NULL,NULL,fem_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
    172175               
    173176                        //decondition" pressure
    … …  
    183186        result=new Result(results->Size()+1,0,1,"p_g",pg);
    184187        results->AddObject(result);
     188
     189        /*output if we have rifts: */
     190        if(numrifts){
     191                OutputRiftsx( &riftproperties,fem_dh->loads,numrifts);
     192                result=new Result(results->Size()+1,0,1,"riftproperties",riftproperties);
     193                results->AddObject(result);
     194        }
    185195}

  • issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp

    r1784 r1805  
    1212#include "./parallel.h"
    1313
    14 void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
     14void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, DataSet* input_loads,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){
    1515
    1616
    … …  
    4747        fem->parameters->FindParam((void*)&solver_string,"solverstring");
    4848
    49         /*Copy loads for backup: */
    50         loads=fem->loads->Copy();
     49        /*Were loads requested as output? : */
     50        if(!input_loads){
     51                loads=fem->loads->Copy(); //we don't want to clobber loads, as they are not needed in output.
     52        }
     53        else{
     54                loads=input_loads; //we are going to modify the loads!
     55        }
    5156
    5257        /*Initialize ug, uf, ug_old and uf_old */
    … …  
    146151        }
    147152
    148         /*Delete loads: */
    149         delete loads;
     153        /*Delete loads only if no ouput was requested: */
     154        if(!input_loads)delete loads;
    150155
    151156        /*clean up*/

  • issm/trunk/src/c/parallel/objectivefunctionC.cpp

    r1209 r1805  
    7777
    7878        //Run diagnostic with updated parameters.
    79         diagnostic_core_nonlinear(&u_g,NULL,NULL,femmodel,inputs,DiagnosticAnalysisEnum(),sub_analysis_type);
     79        diagnostic_core_nonlinear(&u_g,NULL,NULL,NULL,femmodel,inputs,DiagnosticAnalysisEnum(),sub_analysis_type);
    8080        VecToMPISerial(&u_g_double,u_g); VecFree(&u_g);
    8181        inputs->Add("velocity",u_g_double,numberofdofspernode,numberofnodes);

  • issm/trunk/src/c/parallel/parallel.h

    r1671 r1805  
    1818void thermal_core_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
    1919
    20 void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
     20void thermalstatic_core(DataSet* results,FemModel* fems, ParameterInputs* inputs);
     21
     22void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, DataSet* input_loads, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
    2123void diagnostic_core_linear(Vec* ppg,FemModel* fem,ParameterInputs* inputs,int  analysis_type,int sub_analysis_type);
    2224void convergence(int* pconverged, Mat K_ff,Vec p_f,Vec u_f,Vec u_f_old,DataSet* parameters);

  • issm/trunk/src/m/classes/@model/model.m

    r1755 r1805  
    7676        md.rifts=NaN;
    7777        md.riftinfo=NaN;
     78        md.riftproperties=NaN;
    7879        md.numrifts=0;
    7980

  • issm/trunk/src/m/classes/public/loadresultsfromdisk.m

    r1790 r1805  
    1616end
    1717
    18 structure=parseresultsfromdisk(filename);
     18structure=parseresultsfromdisk(filename)
    1919md.results.(structure.analysis_type)=structure;
    2020

  • issm/trunk/src/m/classes/public/plot/plot_manager.m

    r1758 r1805  
    6969                        plot_riftpenetration(md,options_structure,width,i);
    7070                        return;
     71                case 'riftfraction',
     72                        plot_riftfraction(md,options_structure,width,i);
     73                        return;
    7174                case 'sarpwr',
    7275                        plot_sarpwr(md,options_structure,width,i)

  • issm/trunk/src/m/classes/public/plot/plotdoc.m

    r1745 r1805  
    2929disp('                  - ''riftrelvel'': relative velocities along rifts');
    3030disp('                  - ''riftpenetration'': penetration levels for a fault');
     31disp('                  - ''riftfraction'': fill fractions for every node of the rifts');
    3132disp('                  - ''strainrate_tensor'': plot the components of the strainrate tensor (exx,eyy,ezz,exy,exz,eyz) if computed');
    3233disp('                  - ''strainrate_principal'': plot the strainrate tensor principal axis and principal values)');

  • issm/trunk/src/m/classes/public/presolve.m

    r1784 r1805  
    2626        md.riftinfo(count:count+numpairsforthisrift-1,9)=md.rifts(i).friction;
    2727        md.riftinfo(count:count+numpairsforthisrift-1,10)=md.rifts(i).fraction;
     28        md.riftinfo(count:count+numpairsforthisrift-1,11)=md.rifts(i).fractionincrement;
    2829        count=count+numpairsforthisrift;
    2930end

  • issm/trunk/src/m/classes/public/tres.m

    r1787 r1805  
    1313        md.vel=md.results.diagnostic.vel;
    1414        md.pressure=md.results.diagnostic.pressure;
     15        if md.numrifts,
     16                md.riftproperties=md.results.diagnostic.riftproperties;
     17        end
    1518elseif strcmpi(string,'thermalstatic'),
    1619        md.vx=md.results.thermalstatic.vx;

  • issm/trunk/src/m/solutions/cielo/diagnostic_core.m

    r1656 r1805  
    1919ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
    2020isstokes=m_ds.parameters.isstokes;
     21numrifts=m_dhu.parameters.numrifts;
    2122
    2223if ishutter,
    … …  
    5354
    5455        displaystring(debug,'\n%s',['computing horizontal velocities...']);
    55         u_g=diagnostic_core_nonlinear(m_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
     56        [u_g m_dh.loads]=diagnostic_core_nonlinear(m_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
    5657
    5758        displaystring(debug,'\n%s',['computing pressure according to MacAyeal...']);
    … …  
    103104
    104105                displaystring(debug,'\n%s',['computing stokes velocities and pressure ...']);
    105                 u_g=diagnostic_core_nonlinear(m_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
     106                u_g =diagnostic_core_nonlinear(m_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
    106107       
    107108                %"decondition" pressure
    … …  
    115116results.u_g=u_g;
    116117results.p_g=p_g;
     118
     119if numrifts,
     120        results.riftproperties=OutputRifts(m_dh.loads,m_dh.parameters);
     121end

  • issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m

    r1786 r1805  
    8686                varargout(2)={K_fs};
    8787        end
     88        if nout==1,
     89                varargout(1)={loads};
     90        end
    8891end

  • issm/trunk/src/m/solutions/cielo/thermalstatic_core.m

    r1787 r1805  
    4444        results_thermal=thermal_core(models,inputs);
    4545       
    46         %save results from thermal
    47         results.t_g=results_thermal.t_g;
    48         results.m_g=results_thermal.m_g;
    49 
    5046        %add temperature to inputs.
    5147        %Compute depth averaged temperature
    52         temperature_average=FieldDepthAverage(m_t.elements,m_t.nodes,m_t.loads,m_t.materials,results.t_g,'temperature');
     48        temperature_average=FieldDepthAverage(m_t.elements,m_t.nodes,m_t.loads,m_t.materials,results_thermal.t_g,'temperature');
    5349        inputs=add(inputs,'temperature_average',temperature_average,'doublevec',1,m_t.parameters.numberofnodes);
    54         inputs=add(inputs,'temperature',results.t_g,'doublevec',1,m_t.parameters.numberofnodes);
     50        inputs=add(inputs,'temperature',results_thermal.t_g,'doublevec',1,m_t.parameters.numberofnodes);
    5551       
    5652        %now compute diagnostic velocity using the steady state temperature.
    5753        results_diagnostic=diagnostic_core(models,inputs);
    5854       
    59         %save results from thermal
    60         results.u_g=results_diagnostic.u_g;
    61         results.p_g=results_diagnostic.p_g;
    62 
    6355        %convergence?
    64         du_g=results.u_g-u_g_old;
     56        du_g=results_diagnostic.u_g-u_g_old;
    6557        ndu=norm(du_g,2);
    6658        nu=norm(u_g_old,2);
    6759       
    68         dt_g=results.t_g-t_g_old;
     60        dt_g=results_thermal.t_g-t_g_old;
    6961        ndt=norm(dt_g,2);
    7062        nt=norm(t_g_old,2);
    7163
    72         u_g_old=results.u_g;
    73         t_g_old=results.t_g;
     64        u_g_old=results_diagnostic.u_g;
     65        t_g_old=results_thermal.t_g;
    7466       
    7567        displaystring(debug,'%-60s%g\n                                     %s%g\n                                     %s%g%s',...
    … …  
    8779end
    8880
     81%save results from thermal and diagnostic
     82results.u_g=results_diagnostic.u_g;
     83results.p_g=results_diagnostic.p_g;
     84results.t_g=results_thermal.t_g;
     85results.m_g=results_thermal.m_g;
    8986results.step=step;
    9087results.time=0;

  • issm/trunk/src/mex/Makefile.am

    r1665 r1805  
    3030                                NormalizeConstraints\
    3131                                Orth\
     32                                OutputRifts\
    3233                                ParamsToWorkspace\
    3334                                PenaltyConstraints\
    … …  
    157158                          Orth/Orth.h
    158159
     160OutputRifts_SOURCES = OutputRifts/OutputRifts.cpp\
     161                          OutputRifts/OutputRifts.h
     162
    159163ParamsToWorkspace_SOURCES = ParamsToWorkspace/ParamsToWorkspace.cpp\
    160164                          ParamsToWorkspace/ParamsToWorkspace.h

  • issm/trunk/src/mex/TriMeshProcessRifts/TriMeshProcessRifts.cpp

    r1784 r1805  
    2727        /*empty rifts structure: */
    2828        double* pNaN=NULL;
    29         const   char*   fnames[8];
     29        const   char*   fnames[9];
    3030        mwSize     ndim=2;
    3131        mwSize          dimensions[2] = {1,1};
    … …  
    240240                fnames[6] = "friction";
    241241                fnames[7] = "fraction";
     242                fnames[8] = "fractionincrement";
    242243
    243244                dimensions[0]=out_numrifts;
    244245
    245                 pmxa_array=mxCreateStructArray( ndim,dimensions,8,fnames);
     246                pmxa_array=mxCreateStructArray( ndim,dimensions,9,fnames);
    246247               
    247248                for (i=0;i<out_numrifts;i++){
    … …  
    284285                        mxSetField(pmxa_array,i,"penaltypairs",pmxa_array3);
    285286
    286                         /*Friction fraction and fill: set to 0 both */
     287                        /*Friction fraction, fractionincrement  and fill: */
    287288                        mxSetField(pmxa_array,i,"friction",mxCreateDoubleScalar(0));
    288289                        mxSetField(pmxa_array,i,"fill",mxCreateDoubleScalar(IceEnum())); //default is ice
    289290                        mxSetField(pmxa_array,i,"fraction",mxCreateDoubleScalar(0)); //default is ice
     291                        mxSetField(pmxa_array,i,"fractionincrement",mxCreateDoubleScalar(0.1));
    290292                }
    291293        }
Note: See TracChangeset for help on using the changeset viewer.