Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 1804) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 1805) @@ -1383,2 +1383,20 @@ } + +void DataSet::OutputRifts(Vec riftproperties){ + + vector::iterator object; + Riftfront* riftfront=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if((*object)->Enum()==RiftfrontEnum()){ + + riftfront=(Riftfront*)(*object); + riftfront->OutputProperties(riftproperties); + } + } + + + +} Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 1804) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 1805) @@ -82,4 +82,5 @@ void FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse); void UpdateNodePositions(double* thickness,double* bed); + void OutputRifts(Vec riftproperties); }; Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 1804) +++ /issm/trunk/src/c/Makefile.am (revision 1805) @@ -292,5 +292,7 @@ ./FieldDepthAveragex/FieldDepthAveragex.h\ ./FieldExtrudex/FieldExtrudex.cpp\ - ./FieldExtrudex/FieldExtrudex.h + ./FieldExtrudex/FieldExtrudex.h\ + ./OutputRiftsx/OutputRiftsx.h\ + ./OutputRiftsx/OutputRiftsx.cpp @@ -577,4 +579,6 @@ ./FieldExtrudex/FieldExtrudex.cpp\ ./FieldExtrudex/FieldExtrudex.h\ + ./OutputRiftsx/OutputRiftsx.h\ + ./OutputRiftsx/OutputRiftsx.cpp\ ./parallel/diagnostic_core.cpp\ ./parallel/diagnostic_core_linear.cpp\ @@ -596,4 +600,5 @@ ./parallel/transient_core_2d.cpp\ ./parallel/transient_core_3d.cpp\ + ./parallel/thermalstatic_core.cpp\ ./parallel/OutputResults.cpp\ ./parallel/OutputControl.cpp @@ -604,5 +609,5 @@ bin_PROGRAMS = else -bin_PROGRAMS = diagnostic.exe control.exe thermal.exe prognostic.exe transient.exe +bin_PROGRAMS = diagnostic.exe control.exe thermal.exe prognostic.exe transient.exe thermalstatic.exe endif @@ -612,4 +617,7 @@ diagnostic_exe_CXXFLAGS= -fPIC -D_PARALLEL_ +thermalstatic_exe_SOURCES = parallel/thermalstatic.cpp +thermalstatic_exe_CXXFLAGS= -fPIC -D_PARALLEL_ + control_exe_SOURCES = parallel/control.cpp control_exe_CXXFLAGS= -fPIC -D_PARALLEL_ Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 1804) +++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 1805) @@ -18,5 +18,4 @@ void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,Model* model,ConstDataHandle model_handle){ - /*create parameters common to all solutions: */ @@ -125,4 +124,5 @@ throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s"," analysis_type: ",model->analysis_type," sub_analysis_type: ",model->sub_analysis_type," not supported yet!")); } + } Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 1804) +++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 1805) @@ -203,4 +203,12 @@ param->SetInteger(numberofdofspernode); parameters->AddObject(param); + + /*numrifts: */ + ModelFetchData((void**)&model->riftinfo,&model->numrifts,NULL,model_handle,"riftinfo","Matrix","Mat"); + count++; + param= new Param(count,"numrifts",INTEGER); + param->SetInteger(model->numrifts); + parameters->AddObject(param); + xfree((void**)&model->riftinfo); /*All our datasets are already ordered by ids. Set presort flag so that later on, when sorting is requested on these Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 1804) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 1805) @@ -55,4 +55,5 @@ double riftfront_friction; double riftfront_fraction; + double riftfront_fractionincrement; bool riftfront_shelf; double riftfront_penalty_offset; @@ -68,4 +69,5 @@ double friction; double fraction; + double fractionincrement; /*Create loads: */ @@ -186,4 +188,5 @@ friction=*(model->riftinfo+RIFTINFOSIZE*i+8); fraction=*(model->riftinfo+RIFTINFOSIZE*i+9); + fractionincrement=*(model->riftinfo+RIFTINFOSIZE*i+10); strcpy(riftfront_type,"2d"); @@ -209,4 +212,5 @@ riftfront_friction=friction; riftfront_fraction=fraction; + riftfront_fractionincrement=fractionincrement; riftfront_shelf=(bool)model->gridoniceshelf[grid1-1]; @@ -218,5 +222,5 @@ riftfront_prestable=0; - 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); + 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); loads->AddObject(riftfront); Index: /issm/trunk/src/c/ModelProcessorx/Model.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.h (revision 1804) +++ /issm/trunk/src/c/ModelProcessorx/Model.h (revision 1805) @@ -11,5 +11,5 @@ #include "../toolkits/toolkits.h" -#define RIFTINFOSIZE 10 +#define RIFTINFOSIZE 11 struct Model { Index: /issm/trunk/src/c/OutputRiftsx/OutputRiftsx.cpp =================================================================== --- /issm/trunk/src/c/OutputRiftsx/OutputRiftsx.cpp (revision 1805) +++ /issm/trunk/src/c/OutputRiftsx/OutputRiftsx.cpp (revision 1805) @@ -0,0 +1,33 @@ +/*!\file OutputRiftsx + * \brief: output results from diagnostic solution, for rifts. Notably: fraction of + * melange, and penetration. + */ + +#include "./OutputRiftsx.h" + +#undef __FUNCT__ +#define __FUNCT__ "OutputRiftsx" + +#include "../shared/shared.h" +#include "../include/macros.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void OutputRiftsx( Vec* priftproperties, DataSet* loads, int numrifts){ + + /*output: */ + Vec riftproperties=NULL; + + /*Allocate grad_g: */ + riftproperties=NewVec(numrifts); + + /*Compute rift properties : */ + loads->OutputRifts(riftproperties); + + /*Assemble vector: */ + VecAssemblyBegin(riftproperties); + VecAssemblyEnd(riftproperties); + + /*Assign output pointers: */ + *priftproperties=riftproperties; +} Index: /issm/trunk/src/c/OutputRiftsx/OutputRiftsx.h =================================================================== --- /issm/trunk/src/c/OutputRiftsx/OutputRiftsx.h (revision 1805) +++ /issm/trunk/src/c/OutputRiftsx/OutputRiftsx.h (revision 1805) @@ -0,0 +1,14 @@ +/*!\file: OutputRiftsx.h + * \brief header file for rift results output. + */ + +#ifndef _OUTPUTRIFTSX_H +#define _OUTPUTRIFTSX_H + +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void OutputRiftsx( Vec* priftproperties, DataSet* loads, int numrifts); + +#endif /* _OUTPUTRIFTSX_H */ + Index: /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp =================================================================== --- /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp (revision 1804) +++ /issm/trunk/src/c/PenaltyConstraintsx/RiftConstraints.cpp (revision 1805) @@ -16,47 +16,6 @@ int potential; - /*First, has the material non-linearity stabilized? : */ - if(!IsMaterialStable(loads,inputs,analysis_type)){ - /*Do nothing, no constraints activated!: */ - printf("converging material\n"); - num_unstable_constraints=0; - converged=0; - } - else{ - /*Ok, start constraining: */ - printf("converged material: constraining\n"); - - #ifdef _ZIGZAGCOUNTER_ - if(!IsPreStable(loads)){ - - PreConstrain(&num_unstable_constraints,loads,inputs,analysis_type); - - if(!num_unstable_constraints){ - /*Ok, preconstraining is done. Figure out size of unstable constraints pool. Then set - * prestable flag to 1 for all rift grids: */ - - potential=PotentialUnstableConstraints(loads,inputs,analysis_type); - printf(" set of constraints is prestabilised; unstable constraints potentialnumber: %i\n",potential); - - SetPreStable(loads); - num_unstable_constraints=-1; //ensure convergence does not happen! - } - - } - else{ - /*Ok, set of constraints is stable. Start increasing number of inactive constraints, one by one: */ - - /*Figure out max penetration: */ - MaxPenetrationInInputs(loads,inputs,analysis_type); - - /*Constrain rifts: */ - Constrain(&num_unstable_constraints,loads,inputs,analysis_type); - } - #else - Constrain(&num_unstable_constraints,loads,inputs,analysis_type); - #endif - - if(num_unstable_constraints==0)converged=1; - } + Constrain(&num_unstable_constraints,loads,inputs,analysis_type); + if(num_unstable_constraints==0)converged=1; /*Assign output pointers: */ Index: /issm/trunk/src/c/issm.h =================================================================== --- /issm/trunk/src/c/issm.h (revision 1804) +++ /issm/trunk/src/c/issm.h (revision 1805) @@ -57,4 +57,5 @@ #include "./NodeConnectivityx/NodeConnectivityx.h" #include "./ElementConnectivityx/ElementConnectivityx.h" +#include "./OutputRiftsx/OutputRiftsx.h" Index: /issm/trunk/src/c/objects/Riftfront.cpp =================================================================== --- /issm/trunk/src/c/objects/Riftfront.cpp (revision 1804) +++ /issm/trunk/src/c/objects/Riftfront.cpp (revision 1805) @@ -25,5 +25,5 @@ } -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){ +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_fractionincrement, double riftfront_penalty_offset, int riftfront_penalty_lock, bool riftfront_active,int riftfront_counter,bool riftfront_prestable,bool riftfront_shelf){ int i; @@ -54,4 +54,5 @@ friction=riftfront_friction; fraction=riftfront_fraction; + fractionincrement=riftfront_fractionincrement; penalty_offset=riftfront_penalty_offset; penalty_lock=riftfront_penalty_lock; @@ -85,5 +86,5 @@ printf("normal [%g,%g], length %g\n",normal[0],normal[1],normal[2]); - printf("fill: %i friction %g fraction %g\n",fill,friction,fraction); + printf("fill: %i friction %g fraction %g fractionincrement %g \n",fill,friction,fraction,fractionincrement); printf("penalty_offset %g\n",penalty_offset); printf("penalty_lock %i\n",penalty_lock); @@ -114,5 +115,5 @@ printf("normal [%g,%g], length %g\n",normal[0],normal[1],normal[2]); - printf("fill: %i friction %g fraction %g\n",fill,friction,fraction); + printf("fill: %i friction %g fraction %g fractionincrement %g \n",fill,friction,fraction,fractionincrement); printf("penalty_offset %g\n",penalty_offset); printf("penalty_lock %i\n",penalty_lock); @@ -156,4 +157,5 @@ memcpy(marshalled_dataset,&friction,sizeof(friction));marshalled_dataset+=sizeof(friction); memcpy(marshalled_dataset,&fraction,sizeof(fraction));marshalled_dataset+=sizeof(fraction); + memcpy(marshalled_dataset,&fractionincrement,sizeof(fractionincrement));marshalled_dataset+=sizeof(fractionincrement); memcpy(marshalled_dataset,&penalty_offset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset); memcpy(marshalled_dataset,&penalty_lock,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock); @@ -184,4 +186,5 @@ sizeof(friction)+ sizeof(fraction)+ + sizeof(fractionincrement)+ sizeof(penalty_offset)+ sizeof(penalty_lock)+ @@ -225,4 +228,5 @@ memcpy(&friction,marshalled_dataset,sizeof(friction));marshalled_dataset+=sizeof(friction); memcpy(&fraction,marshalled_dataset,sizeof(fraction));marshalled_dataset+=sizeof(fraction); + memcpy(&fractionincrement,marshalled_dataset,sizeof(fractionincrement));marshalled_dataset+=sizeof(fractionincrement); memcpy(&penalty_offset,marshalled_dataset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset); memcpy(&penalty_lock,marshalled_dataset,sizeof(penalty_lock));marshalled_dataset+=sizeof(penalty_lock); @@ -621,49 +625,18 @@ penetration=(vxvy_list[1][0]-vxvy_list[0][0])*normal[0]+(vxvy_list[1][1]-vxvy_list[0][1])*normal[1]; - #ifdef _ZIGZAGCOUNTER_ - - found=inputs->Recover("max_penetration",&max_penetration); - if(!found)throw ErrorException(__FUNCT__," could not find max_penetration in inputs!"); - - /*Activate or deactivate penalties: */ - if(penetration<0){ - //printf("riftfront %i penetrating.\n",this->GetId()); - /*There is penetration, we need to active the penalty so this penetration will be NULL: */ - activate=1; - } - else{ - /*Ok, there is no penetration for these two grids of the rift. We want to deactivate the penalty. If we do - * it all at once, then we will zigzag forever. Only deactivate once at a time. Deactivate the one that most - * wants to, ie the one with the max penetration: */ - if (penetration==max_penetration){ - //printf("riftfront %i at max penetration.\n",this->GetId()); - activate=0; - } - else{ - /*Ok, we are dealing with the rest of the penalties that want to be freed. But may be in this lot there - * are penalties that were already free? Keep them as such. For the rest, do not release them yet: */ - if (this->active==0){ - activate=0; - } - else{ - activate=1; - } - } - } - - /*Here, we try to avoid zigzaging. When a penalty activates and deactivates for more than penalty_lock times, - * we lock it: */ - if(this->counter>this->penalty_lock){ - /*This penalty has zig zagged too long, fix it to smooth results: */ - activate=1; this->active=activate; - if (this->counter==(this->penalty_lock+1)){ - printf("riftfront %i locking at %g m/yr\n",this->GetId(),penetration*365.0*24.0*3600.0); - this->counter++; - } - } - #else + /*activation: */ if(penetration<0)activate=1; else activate=0; - #endif + + /*Here, we try to avoid zigzaging. When a penalty activates and deactivates for more than penalty_lock times, + * we increase the fraction of melange:*/ + if(this->counter>this->penalty_lock){ + /*reset counter: */ + this->counter=0; + /*increase melange fraction: */ + this->fraction+=this->fractionincrement; + if (this->fraction>1)this->fraction=(double)1.0; + //printf("riftfront %i fraction: %g\n",this->GetId(),this->fraction); + } //Figure out stability of this penalty @@ -814,2 +787,17 @@ return this->material_converged; } + +#undef __FUNCT__ +#define __FUNCT__ "Riftfront::OutputProperties" +void Riftfront::OutputProperties(Vec riftproperties){ + + int row_id=0; + double value; + + /*recover id of penalty: */ + row_id=this->GetId()-1; //c indexing, ids were matlab indexed + value=(double)this->fraction; + + /*Plug id and fraction into riftproperties matrix: */ + VecSetValues(riftproperties,1,&row_id,&value,INSERT_VALUES); +} Index: /issm/trunk/src/c/objects/Riftfront.h =================================================================== --- /issm/trunk/src/c/objects/Riftfront.h (revision 1804) +++ /issm/trunk/src/c/objects/Riftfront.h (revision 1805) @@ -42,4 +42,5 @@ double friction; double fraction; + double fractionincrement; bool shelf; @@ -57,5 +58,5 @@ Riftfront(); - 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); + 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); ~Riftfront(); @@ -84,4 +85,5 @@ int PotentialUnstableConstraint(int* punstable, void* inputs, int analysis_type); int IsMaterialStable(void* inputs, int analysis_type); + void OutputProperties(Vec riftproperties); Object* copy(); }; Index: /issm/trunk/src/c/parallel/ControlInitialization.cpp =================================================================== --- /issm/trunk/src/c/parallel/ControlInitialization.cpp (revision 1804) +++ /issm/trunk/src/c/parallel/ControlInitialization.cpp (revision 1805) @@ -97,5 +97,5 @@ //horizontal velocity if(debug)_printf_("%s\n"," computing horizontal velocities..."); - diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); + diagnostic_core_nonlinear(&ug,NULL,NULL,NULL, fem_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); if(debug)_printf_("%s\n"," extruding horizontal velocities..."); VecDuplicatePatch(&ug_horiz,ug); FieldExtrudex( ug_horiz,fem_dh->elements,fem_dh->nodes, fem_dh->loads,fem_dh-> materials,"velocity",1); @@ -138,5 +138,5 @@ if(debug)_printf_("%s\n"," computing stokes velocities and pressure ..."); VecFree(&ug); - diagnostic_core_nonlinear(&ug,NULL,NULL,fem_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); + diagnostic_core_nonlinear(&ug,NULL,NULL,NULL, fem_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); VecFree(&pg); Index: /issm/trunk/src/c/parallel/GradJCompute.cpp =================================================================== --- /issm/trunk/src/c/parallel/GradJCompute.cpp (revision 1804) +++ /issm/trunk/src/c/parallel/GradJCompute.cpp (revision 1805) @@ -59,5 +59,5 @@ _printf_("%s\n"," recover solution for this stiffness and right hand side:"); - diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,femmodel,inputs,DiagnosticAnalysisEnum(),sub_analysis_type); + diagnostic_core_nonlinear(&u_g,&K_ff0,&K_fs0,NULL, femmodel,inputs,DiagnosticAnalysisEnum(),sub_analysis_type); VecToMPISerial(&u_g_double,u_g); VecFree(&u_g); inputs->Add("velocity",u_g_double,numberofdofspernode,numberofnodes); Index: /issm/trunk/src/c/parallel/ProcessResults.cpp =================================================================== --- /issm/trunk/src/c/parallel/ProcessResults.cpp (revision 1804) +++ /issm/trunk/src/c/parallel/ProcessResults.cpp (revision 1805) @@ -70,4 +70,7 @@ double* param_g=NULL; double* parameter=NULL; + Vec riftproperties=NULL; + double* riftproperties_serial=NULL; + int numrifts=0; Vec s_g=NULL; @@ -106,4 +109,18 @@ fem_dhu=fems+3; fem_sl=fems+4; + + /*some flags needed: */ + fem_dh->parameters->FindParam((void*)&dim,"dim"); + fem_dhu->parameters->FindParam((void*)&ishutter,"ishutter"); + fem_ds->parameters->FindParam((void*)&isstokes,"isstokes"); + fem_dh->parameters->FindParam((void*)&ismacayealpattyn,"ismacayealpattyn"); + } + if(analysis_type==ThermalstaticAnalysisEnum()){ + fem_dh=fems+0; + fem_dv=fems+1; + fem_ds=fems+2; + fem_dhu=fems+3; + fem_sl=fems+4; + fem_t=fems+5; /*some flags needed: */ @@ -411,4 +428,14 @@ xfree((void**)&partition); } + else if(strcmp(result->GetFieldName(),"riftproperties")==0){ + result->GetField(&riftproperties); + fem_dh->parameters->FindParam((void*)&numrifts,"numrifts"); + VecToMPISerial(&riftproperties_serial,riftproperties); + + /*Ok, add parameter to newresults: */ + newresult=new Result(newresults->Size()+1,result->GetTime(),result->GetStep(),"riftproperties",riftproperties_serial,numrifts); + newresults->AddObject(newresult); + + } else{ /*Just copy the result into the new results dataset: */ Index: /issm/trunk/src/c/parallel/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core.cpp (revision 1804) +++ /issm/trunk/src/c/parallel/diagnostic_core.cpp (revision 1805) @@ -35,4 +35,5 @@ Vec slopex=NULL; Vec slopey=NULL; + Vec riftproperties=NULL; /*flags: */ @@ -46,4 +47,5 @@ int numberofdofspernode_ds; int numberofnodes; + int numrifts=0; double stokesreconditioning; @@ -70,4 +72,5 @@ fem_ds->parameters->FindParam((void*)&isstokes,"isstokes"); fem_ds->parameters->FindParam((void*)&stokesreconditioning,"stokesreconditioning"); + fem_dh->parameters->FindParam((void*)&numrifts,"numrifts"); //specific parameters for specific models @@ -112,5 +115,5 @@ if(debug)_printf_("%s\n"," computing horizontal velocities..."); - diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); + diagnostic_core_nonlinear(&ug,NULL,NULL,fem_dh->loads,fem_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); if(debug)_printf_("%s\n"," computing pressure according to MacAyeal..."); @@ -169,5 +172,5 @@ if(debug)_printf_("%s\n"," computing stokes velocities and pressure ..."); VecFree(&ug); - diagnostic_core_nonlinear(&ug,NULL,NULL,fem_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); + diagnostic_core_nonlinear(&ug,NULL,NULL,NULL,fem_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); //decondition" pressure @@ -183,3 +186,10 @@ result=new Result(results->Size()+1,0,1,"p_g",pg); results->AddObject(result); + + /*output if we have rifts: */ + if(numrifts){ + OutputRiftsx( &riftproperties,fem_dh->loads,numrifts); + result=new Result(results->Size()+1,0,1,"riftproperties",riftproperties); + results->AddObject(result); + } } Index: /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 1804) +++ /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 1805) @@ -12,5 +12,5 @@ #include "./parallel.h" -void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ +void diagnostic_core_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, DataSet* input_loads,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ @@ -47,6 +47,11 @@ fem->parameters->FindParam((void*)&solver_string,"solverstring"); - /*Copy loads for backup: */ - loads=fem->loads->Copy(); + /*Were loads requested as output? : */ + if(!input_loads){ + loads=fem->loads->Copy(); //we don't want to clobber loads, as they are not needed in output. + } + else{ + loads=input_loads; //we are going to modify the loads! + } /*Initialize ug, uf, ug_old and uf_old */ @@ -146,6 +151,6 @@ } - /*Delete loads: */ - delete loads; + /*Delete loads only if no ouput was requested: */ + if(!input_loads)delete loads; /*clean up*/ Index: /issm/trunk/src/c/parallel/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 1804) +++ /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 1805) @@ -77,5 +77,5 @@ //Run diagnostic with updated parameters. - diagnostic_core_nonlinear(&u_g,NULL,NULL,femmodel,inputs,DiagnosticAnalysisEnum(),sub_analysis_type); + diagnostic_core_nonlinear(&u_g,NULL,NULL,NULL,femmodel,inputs,DiagnosticAnalysisEnum(),sub_analysis_type); VecToMPISerial(&u_g_double,u_g); VecFree(&u_g); inputs->Add("velocity",u_g_double,numberofdofspernode,numberofnodes); Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 1804) +++ /issm/trunk/src/c/parallel/parallel.h (revision 1805) @@ -18,5 +18,7 @@ void thermal_core_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); -void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); +void thermalstatic_core(DataSet* results,FemModel* fems, ParameterInputs* inputs); + +void 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); void diagnostic_core_linear(Vec* ppg,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); void convergence(int* pconverged, Mat K_ff,Vec p_f,Vec u_f,Vec u_f_old,DataSet* parameters); Index: /issm/trunk/src/c/parallel/thermalstatic.cpp =================================================================== --- /issm/trunk/src/c/parallel/thermalstatic.cpp (revision 1805) +++ /issm/trunk/src/c/parallel/thermalstatic.cpp (revision 1805) @@ -0,0 +1,144 @@ +/*!\file: thermalstatic.cpp + * \brief: thermalstatic solution + */ + +#include "../issm.h" +#include "./parallel.h" + +#undef __FUNCT__ +#define __FUNCT__ "thermalstatic" + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + int numberofnodes; + int qmu_analysis=0; + + /*Fem models : */ + FemModel femmodels[7]; + + /*Results: */ + DataSet* results=NULL; + Result* result=NULL; + + ParameterInputs* inputs=NULL; + int waitonlock=0; + + double* u_g_initial=NULL; + double* p_g_initial=NULL; + double dt; + Param* param=NULL; + + MODULEBOOT(); + + #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_)) + throw ErrorException(__FUNCT__," parallel executable was compiled without support of parallel libraries!"); + #endif + + PetscInitialize(&argc,&argv,(char *)0,""); + + /*Size and rank: */ + MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); + MPI_Comm_size(MPI_COMM_WORLD,&num_procs); + + _printf_("recover , input file name and output file name:\n"); + inputfilename=argv[2]; + outputfilename=argv[3]; + lockname=argv[4]; + + /*Open handle to data on disk: */ + fid=pfopen(inputfilename,"rb"); + + _printf_("read and create finite element model:\n"); + _printf_("\n reading diagnostic horiz model data:\n"); + CreateFemModel(&femmodels[0],fid,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); + _printf_("\n reading diagnostic vert model data:\n"); + CreateFemModel(&femmodels[1],fid,DiagnosticAnalysisEnum(),VertAnalysisEnum()); + _printf_("\n reading diagnostic stokes model data:\n"); + CreateFemModel(&femmodels[2],fid,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); + _printf_("\n reading diagnostic hutter model data:\n"); + CreateFemModel(&femmodels[3],fid,DiagnosticAnalysisEnum(),HutterAnalysisEnum()); + _printf_("\n reading surface and bed slope computation model data:\n"); + CreateFemModel(&femmodels[4],fid,SlopeComputeAnalysisEnum(),NoneAnalysisEnum()); + _printf_("\n read and create thermal finite element model:\n"); + CreateFemModel(&femmodels[5],fid,ThermalAnalysisEnum(),SteadyAnalysisEnum()); + _printf_("\n read and create melting finite element model:\n"); + CreateFemModel(&femmodels[6],fid,MeltingAnalysisEnum(),SteadyAnalysisEnum()); + + _printf_("initialize inputs:\n"); + femmodels[0].parameters->FindParam((void*)&u_g_initial,"u_g"); + femmodels[0].parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + femmodels[5].parameters->FindParam((void*)&dt,"dt"); + femmodels[5].parameters->FindParam((void*)&p_g_initial,"p_g"); + + inputs=new ParameterInputs; + inputs->Add("velocity",u_g_initial,3,numberofnodes); + inputs->Add("pressure",p_g_initial,1,numberofnodes); + inputs->Add("dt",dt); + + /*erase velocities: */ + param=(Param*)femmodels[0].parameters->FindParamObject("u_g"); + femmodels[0].parameters->DeleteObject((Object*)param); + + _printf_("initialize results:\n"); + results=new DataSet(ResultsEnum()); + + /*are we running the solution sequence, or a qmu wrapper around it? : */ + femmodels[0].parameters->FindParam((void*)&qmu_analysis,"qmu_analysis"); + if(!qmu_analysis){ + + /*run diagnostic analysis: */ + _printf_("call computational core:\n"); + thermalstatic_core(results,femmodels,inputs); + + } + else{ + /*run qmu analysis: */ + _printf_("calling qmu analysis on thermalstatic core:\n"); + + #ifdef _HAVE_DAKOTA_ + Qmux(&femmodels[0],inputs,ThermalstaticAnalysisEnum(),NoneAnalysisEnum()); + #else + throw ErrorException(__FUNCT__," Dakota not present, cannot do qmu!"); + #endif + } + + /*Add analysis_type to results: */ + result=new Result(results->Size()+1,0,1,"analysis_type","thermalstatic"); + results->AddObject(result); + + _printf_("process results:\n"); + ProcessResults(&results,&femmodels[0],ThermalstaticAnalysisEnum()); + + _printf_("write results to disk:\n"); + OutputResults(results,outputfilename); + + _printf_("write lock file:\n"); + femmodels[0].parameters->FindParam((void*)&waitonlock,"waitonlock"); + if (waitonlock){ + WriteLockFile(lockname); + } + _printf_("closing MPI and Petsc\n"); + PetscFinalize(); + + + /*end module: */ + MODULEEND(); + + /*Free ressources */ + xfree((void**)&u_g_initial); + xfree((void**)&p_g_initial); + + return 0; //unix success return; +} Index: /issm/trunk/src/c/parallel/thermalstatic_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/thermalstatic_core.cpp (revision 1805) +++ /issm/trunk/src/c/parallel/thermalstatic_core.cpp (revision 1805) @@ -0,0 +1,148 @@ +/*!\file: thermalstatic_core.cpp + * \brief: core of the thermalstatic solution + */ + +#undef __FUNCT__ +#define __FUNCT__ "thermalstatic_core" + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../issm.h" + +void thermalstatic_core(DataSet* results,FemModel* fems, ParameterInputs* inputs){ + + extern int my_rank; + + /*fem models: */ + FemModel* fem_dh=NULL; + FemModel* fem_dv=NULL; + FemModel* fem_dhu=NULL; + FemModel* fem_ds=NULL; + FemModel* fem_sl=NULL; + FemModel* fem_t=NULL; + FemModel* fem_m=NULL; + + /*output: */ + Result* result=NULL; + DataSet* results_thermal=NULL; + DataSet* results_diagnostic=NULL; + + /*solutions: */ + Vec u_g=NULL; + Vec old_u_g=NULL; + Vec t_g=NULL; + Vec t_g_average=NULL; + Vec old_t_g=NULL; + Vec p_g=NULL; + Vec m_g=NULL; + Vec du_g=NULL; + Vec dt_g=NULL; + double ndu,nu; + double ndt,nt; + double eps_rel; + + /*flags: */ + int debug=0; + int isstokes=0; + int numberofnodes; + int ndof; + int converged; + int step; + + /*recover fem models: */ + fem_dh=fems+0; + fem_dv=fems+1; + fem_ds=fems+2; + fem_dhu=fems+3; + fem_sl=fems+4; + fem_t=fems+5; + fem_m=fems+6; + + //first recover parameters common to all solutions + fem_dh->parameters->FindParam((void*)&debug,"debug");debug=1; + fem_dh->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + fem_dh->parameters->FindParam((void*)&eps_rel,"eps_rel"); + fem_ds->parameters->FindParam((void*)&isstokes,"isstokes"); + + + + //initialize: + converged=0; + step=1; + + if (isstokes)ndof=4; + else ndof=3; + + for(;;){ + + if(debug)_printf_("%s%i\n","computing temperature and velocity for step: ",step); + + //first compute temperature at steady state. + if (step>1){ + inputs->Add("velocity",u_g,ndof,numberofnodes); + } + results_thermal=new DataSet(ResultsEnum()); + thermal_core(results_thermal,fems+5,inputs); + + //get t_g and m_g; + VecFree(&t_g);results_thermal->FindResult(&t_g,"t_g"); + VecFree(&m_g);results_thermal->FindResult(&m_g,"m_g"); + delete results_thermal; + + //Add temperature to inputs. + //compute depth averaged temperature and add to inputs + VecDuplicatePatch(&t_g_average,t_g); + FieldDepthAveragex( t_g_average, fem_t->elements,fem_t->nodes, fem_t->loads, fem_t->materials,"temperature"); + inputs->Add("temperature_average",t_g_average,1,numberofnodes); + inputs->Add("temperature",t_g,1,numberofnodes); + VecFree(&t_g_average); //not needed anymore + + //now compute diagnostic velocity using the steady state temperature. + results_diagnostic=new DataSet(ResultsEnum()); + diagnostic_core(results_diagnostic,fems, inputs); + + //get p_g and u_g + VecFree(&u_g);results_diagnostic->FindResult(&u_g,"u_g"); + VecFree(&p_g);results_diagnostic->FindResult(&p_g,"p_g"); + delete results_diagnostic; + + //convergence? + if(step>1){ + VecFree(&du_g);VecDuplicatePatch(&du_g,old_u_g);VecAYPX(du_g,-1.0,u_g); + VecNorm(du_g,NORM_2,&ndu); VecNorm(old_u_g,NORM_2,&nu); + + VecFree(&dt_g);VecDuplicatePatch(&dt_g,old_t_g); VecAYPX(dt_g,-1.0,t_g); + VecNorm(dt_g,NORM_2,&ndu); VecNorm(old_t_g,NORM_2,&nu); + + + if (debug) _printf_("%-60s%g\n %s%g\n %s%g%s\n", + " relative convergence criterion: velocity -> norm(du)/norm(u)= ",ndu/nu*100," temperature -> norm(dt)/norm(t)=",ndt/nt*100," eps_rel: ",eps_rel*100," %"); + + if ((ndu/nu<=eps_rel) && (ndt/nt<=eps_rel)) converged=1; + else converged=0; + } + else{ + converged=0; + } + + VecFree(&old_u_g);VecDuplicatePatch(&old_u_g,u_g); + VecFree(&old_t_g);VecDuplicatePatch(&old_t_g,t_g); + + step++; + if (converged)break; + } + + + /*Plug results into output dataset: */ + result=new Result(results->Size()+1,0,1,"u_g",u_g); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"p_g",p_g); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"t_g",t_g); + results->AddObject(result); + result=new Result(results->Size()+1,0,1,"m_g",m_g); + results->AddObject(result); +} Index: /issm/trunk/src/m/classes/@model/model.m =================================================================== --- /issm/trunk/src/m/classes/@model/model.m (revision 1804) +++ /issm/trunk/src/m/classes/@model/model.m (revision 1805) @@ -76,4 +76,5 @@ md.rifts=NaN; md.riftinfo=NaN; + md.riftproperties=NaN; md.numrifts=0; Index: /issm/trunk/src/m/classes/public/loadresultsfromdisk.m =================================================================== --- /issm/trunk/src/m/classes/public/loadresultsfromdisk.m (revision 1804) +++ /issm/trunk/src/m/classes/public/loadresultsfromdisk.m (revision 1805) @@ -16,5 +16,5 @@ end -structure=parseresultsfromdisk(filename); +structure=parseresultsfromdisk(filename) md.results.(structure.analysis_type)=structure; Index: /issm/trunk/src/m/classes/public/plot/plot_manager.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_manager.m (revision 1804) +++ /issm/trunk/src/m/classes/public/plot/plot_manager.m (revision 1805) @@ -69,4 +69,7 @@ plot_riftpenetration(md,options_structure,width,i); return; + case 'riftfraction', + plot_riftfraction(md,options_structure,width,i); + return; case 'sarpwr', plot_sarpwr(md,options_structure,width,i) Index: /issm/trunk/src/m/classes/public/plot/plot_riftfraction.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plot_riftfraction.m (revision 1805) +++ /issm/trunk/src/m/classes/public/plot/plot_riftfraction.m (revision 1805) @@ -0,0 +1,53 @@ +function plot_riftfraction(md,options_structure,width,i); +%PLOT_RIFTRELVEL - plot rift fractions +% +% Usage: +% plot_riftfraction(md,options_structure,width,i); +% +% See also: PLOTMODEL + +%plot mesh boundaries +subplot(width,width,i); + +%units +if ~isnan(options_structure.unitmultiplier), + md.x=md.x*options_structure.unitmultiplier; + md.y=md.y*options_structure.unitmultiplier; + md.z=md.z*options_structure.unitmultiplier; +end + +for i=1:size(md.segments,1), + h1=plot(md.x(md.segments(i,1:2)),md.y(md.segments(i,1:2)),'k.-');hold on; +end + +%first, build a vector of fractions, over all grids. +fractions=zeros(md.numberofgrids,1); + +%plug riftproperties fractions: +fractions(md.riftinfo(:,1))=md.riftproperties; +fractions(md.riftinfo(:,2))=md.riftproperties; + +%complete the tips. +for i=1:length(md.rifts), + tips=md.rifts(i).tips; + fractions(tips)=1; +end + +hold on; +for i=1:length(md.rifts), + segments=md.rifts(i).segments(:,1:2)'; + xc=md.x(segments(:)); + yc=md.y(segments(:)); + zc=fractions(segments(:)); + h2=patch('Xdata',xc,'Ydata',yc,'Zdata',zc,'Cdata',zc,'facecolor','none','edgecolor','flat'); +end +legend([h1,h2],'mesh boundaries','rifts') + +%apply options +if isnan(options_structure.title) + options_structure.title='Rift relative velocities'; +end +if isnan(options_structure.colorbar) + options_structure.colorbar=0; +end +applyoptions(md,[],options_structure); Index: /issm/trunk/src/m/classes/public/plot/plotdoc.m =================================================================== --- /issm/trunk/src/m/classes/public/plot/plotdoc.m (revision 1804) +++ /issm/trunk/src/m/classes/public/plot/plotdoc.m (revision 1805) @@ -29,4 +29,5 @@ disp(' - ''riftrelvel'': relative velocities along rifts'); disp(' - ''riftpenetration'': penetration levels for a fault'); +disp(' - ''riftfraction'': fill fractions for every node of the rifts'); disp(' - ''strainrate_tensor'': plot the components of the strainrate tensor (exx,eyy,ezz,exy,exz,eyz) if computed'); disp(' - ''strainrate_principal'': plot the strainrate tensor principal axis and principal values)'); Index: /issm/trunk/src/m/classes/public/presolve.m =================================================================== --- /issm/trunk/src/m/classes/public/presolve.m (revision 1804) +++ /issm/trunk/src/m/classes/public/presolve.m (revision 1805) @@ -26,4 +26,5 @@ md.riftinfo(count:count+numpairsforthisrift-1,9)=md.rifts(i).friction; md.riftinfo(count:count+numpairsforthisrift-1,10)=md.rifts(i).fraction; + md.riftinfo(count:count+numpairsforthisrift-1,11)=md.rifts(i).fractionincrement; count=count+numpairsforthisrift; end Index: /issm/trunk/src/m/classes/public/tres.m =================================================================== --- /issm/trunk/src/m/classes/public/tres.m (revision 1804) +++ /issm/trunk/src/m/classes/public/tres.m (revision 1805) @@ -13,4 +13,7 @@ md.vel=md.results.diagnostic.vel; md.pressure=md.results.diagnostic.pressure; + if md.numrifts, + md.riftproperties=md.results.diagnostic.riftproperties; + end elseif strcmpi(string,'thermalstatic'), md.vx=md.results.thermalstatic.vx; Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/diagnostic_core.m (revision 1804) +++ /issm/trunk/src/m/solutions/cielo/diagnostic_core.m (revision 1805) @@ -19,4 +19,5 @@ ismacayealpattyn=m_dh.parameters.ismacayealpattyn; isstokes=m_ds.parameters.isstokes; +numrifts=m_dhu.parameters.numrifts; if ishutter, @@ -53,5 +54,5 @@ displaystring(debug,'\n%s',['computing horizontal velocities...']); - u_g=diagnostic_core_nonlinear(m_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); + [u_g m_dh.loads]=diagnostic_core_nonlinear(m_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); displaystring(debug,'\n%s',['computing pressure according to MacAyeal...']); @@ -103,5 +104,5 @@ displaystring(debug,'\n%s',['computing stokes velocities and pressure ...']); - u_g=diagnostic_core_nonlinear(m_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); + u_g =diagnostic_core_nonlinear(m_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); %"decondition" pressure @@ -115,2 +116,6 @@ results.u_g=u_g; results.p_g=p_g; + +if numrifts, + results.riftproperties=OutputRifts(m_dh.loads,m_dh.parameters); +end Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m (revision 1804) +++ /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m (revision 1805) @@ -86,3 +86,6 @@ varargout(2)={K_fs}; end + if nout==1, + varargout(1)={loads}; + end end Index: /issm/trunk/src/m/solutions/cielo/thermalstatic_core.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/thermalstatic_core.m (revision 1804) +++ /issm/trunk/src/m/solutions/cielo/thermalstatic_core.m (revision 1805) @@ -44,32 +44,24 @@ results_thermal=thermal_core(models,inputs); - %save results from thermal - results.t_g=results_thermal.t_g; - results.m_g=results_thermal.m_g; - %add temperature to inputs. %Compute depth averaged temperature - temperature_average=FieldDepthAverage(m_t.elements,m_t.nodes,m_t.loads,m_t.materials,results.t_g,'temperature'); + temperature_average=FieldDepthAverage(m_t.elements,m_t.nodes,m_t.loads,m_t.materials,results_thermal.t_g,'temperature'); inputs=add(inputs,'temperature_average',temperature_average,'doublevec',1,m_t.parameters.numberofnodes); - inputs=add(inputs,'temperature',results.t_g,'doublevec',1,m_t.parameters.numberofnodes); + inputs=add(inputs,'temperature',results_thermal.t_g,'doublevec',1,m_t.parameters.numberofnodes); %now compute diagnostic velocity using the steady state temperature. results_diagnostic=diagnostic_core(models,inputs); - %save results from thermal - results.u_g=results_diagnostic.u_g; - results.p_g=results_diagnostic.p_g; - %convergence? - du_g=results.u_g-u_g_old; + du_g=results_diagnostic.u_g-u_g_old; ndu=norm(du_g,2); nu=norm(u_g_old,2); - dt_g=results.t_g-t_g_old; + dt_g=results_thermal.t_g-t_g_old; ndt=norm(dt_g,2); nt=norm(t_g_old,2); - u_g_old=results.u_g; - t_g_old=results.t_g; + u_g_old=results_diagnostic.u_g; + t_g_old=results_thermal.t_g; displaystring(debug,'%-60s%g\n %s%g\n %s%g%s',... @@ -87,4 +79,9 @@ end +%save results from thermal and diagnostic +results.u_g=results_diagnostic.u_g; +results.p_g=results_diagnostic.p_g; +results.t_g=results_thermal.t_g; +results.m_g=results_thermal.m_g; results.step=step; results.time=0; Index: /issm/trunk/src/mex/Makefile.am =================================================================== --- /issm/trunk/src/mex/Makefile.am (revision 1804) +++ /issm/trunk/src/mex/Makefile.am (revision 1805) @@ -30,4 +30,5 @@ NormalizeConstraints\ Orth\ + OutputRifts\ ParamsToWorkspace\ PenaltyConstraints\ @@ -157,4 +158,7 @@ Orth/Orth.h +OutputRifts_SOURCES = OutputRifts/OutputRifts.cpp\ + OutputRifts/OutputRifts.h + ParamsToWorkspace_SOURCES = ParamsToWorkspace/ParamsToWorkspace.cpp\ ParamsToWorkspace/ParamsToWorkspace.h Index: /issm/trunk/src/mex/OutputRifts/OutputRifts.cpp =================================================================== --- /issm/trunk/src/mex/OutputRifts/OutputRifts.cpp (revision 1805) +++ /issm/trunk/src/mex/OutputRifts/OutputRifts.cpp (revision 1805) @@ -0,0 +1,49 @@ +/*\file OutputRifts.c + *\brief: output rift properties (fraction, penetration, etc ...) + */ + +#include "./OutputRifts.h" + +void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){ + + /*diverse: */ + int noerr=1; + int numrifts; + + /*input datasets: */ + DataSet* loads=NULL; + + /* output datasets: */ + Vec riftproperties=NULL; + + + /*Boot module: */ + MODULEBOOT(); + + /*checks on arguments on the matlab side: */ + CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&OutputRiftsUsage); + + /*Input datasets: */ + FetchData((void**)&loads,NULL,NULL,LOADS,"DataSet",NULL); + FetchData((void**)&numrifts,NULL,NULL,mxGetField(PARAMETERS,0,"numrifts"),"Integer",NULL); + + /*!Call core code: */ + OutputRiftsx(&riftproperties,loads,numrifts); + + /*write output : */ + WriteData(RIFTPROPERTIES,riftproperties,0,0,"Vector",NULL); + + /*Free ressources: */ + delete loads; + VecFree(&riftproperties); + + /*end module: */ + MODULEEND(); +} + +void OutputRiftsUsage(void) +{ + _printf_("\n"); + _printf_(" usage: [riftproperties] = %s(loads,parameters);\n",__FUNCT__); + _printf_("\n"); +} Index: /issm/trunk/src/mex/OutputRifts/OutputRifts.h =================================================================== --- /issm/trunk/src/mex/OutputRifts/OutputRifts.h (revision 1805) +++ /issm/trunk/src/mex/OutputRifts/OutputRifts.h (revision 1805) @@ -0,0 +1,35 @@ + +/* + OutputRifts.h +*/ + + +#ifndef _OUTPUTRIFTS_H +#define _OUTPUTRIFTS_H + +/* local prototypes: */ +void OutputRiftsUsage(void); + +#include "../../c/issm.h" + +#undef __FUNCT__ +#define __FUNCT__ "OutputRifts" + +/* serial input macros: */ +#define LOADS (mxArray*)prhs[0] +#define PARAMETERS (mxArray*)prhs[1] + +/* serial output macros: */ +#define RIFTPROPERTIES (mxArray**)&plhs[0] + +/* serial arg counts: */ +#undef NLHS +#define NLHS 1 +#undef NRHS +#define NRHS 2 + + +#endif /* _OUTPUTRIFTS_H */ + + + Index: /issm/trunk/src/mex/TriMeshProcessRifts/TriMeshProcessRifts.cpp =================================================================== --- /issm/trunk/src/mex/TriMeshProcessRifts/TriMeshProcessRifts.cpp (revision 1804) +++ /issm/trunk/src/mex/TriMeshProcessRifts/TriMeshProcessRifts.cpp (revision 1805) @@ -27,5 +27,5 @@ /*empty rifts structure: */ double* pNaN=NULL; - const char* fnames[8]; + const char* fnames[9]; mwSize ndim=2; mwSize dimensions[2] = {1,1}; @@ -240,8 +240,9 @@ fnames[6] = "friction"; fnames[7] = "fraction"; + fnames[8] = "fractionincrement"; dimensions[0]=out_numrifts; - pmxa_array=mxCreateStructArray( ndim,dimensions,8,fnames); + pmxa_array=mxCreateStructArray( ndim,dimensions,9,fnames); for (i=0;i