Changeset 2110


Ignore:
Timestamp:
09/04/09 11:33:42 (15 years ago)
Author:
Eric.Larour
Message:

Added MassFlux capabilities.

Location:
issm/trunk
Files:
14 added
30 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/config.h.in

    r1385 r2110  
    3939/* Define to 1 if the X Window System is missing or not being used. */
    4040#undef X_DISPLAY_MISSING
    41 
    42 /* Macro to enable debugging in Dakota. */
    43 #undef _DEBUG_
    4441
    4542/* with Blacs in ISSM src */
    … …  
    9087/* with Triangle in ISSM src */
    9188#undef _HAVE_TRIANGLE_
     89
     90/* Macro to enable debugging in ISSM. */
     91#undef _ISSM_DEBUG_

  • issm/trunk/etc/cluster.rc

    r2042 r2110  
    4343
    4444cluster_name=larsen
    45 cluster_codepath=/u/astrid1/larour/issm/trunk/bin
     45cluster_codepath=/u/astrid-r1b/larour/issm/trunk/bin
    4646cluster_executionpath=/u/wilkes-r1b/larour/Testing/Execution
    4747

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

    r2047 r2110  
    274274                                        ./ElementConnectivityx/ElementConnectivityx.cpp\
    275275                                        ./ElementConnectivityx/ElementConnectivityx.h\
     276                                        ./MassFluxx/MassFluxx.cpp\
     277                                        ./MassFluxx/MassFluxx.h\
    276278                                        ./SystemMatricesx/SystemMatricesx.cpp\
    277279                                        ./SystemMatricesx/SystemMatricesx.h\
    … …  
    562564                                        ./NormalizeConstraintsx/NormalizeConstraintsx.cpp\
    563565                                        ./NormalizeConstraintsx/NormalizeConstraintsx.h\
     566                                        ./MassFluxx/MassFluxx.cpp\
     567                                        ./MassFluxx/MassFluxx.h\
    564568                                        ./SystemMatricesx/SystemMatricesx.cpp\
    565569                                        ./SystemMatricesx/SystemMatricesx.h\

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

    r1907 r2110  
    1919void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
    2020
    21 
    2221        /*create parameters common to all solutions: */
    2322        CreateParameters(pparameters,iomodel,iomodel_handle);
    24         CreateParametersQmu(pparameters,iomodel,iomodel_handle);
    2523        CreateParametersControl(pparameters,iomodel,iomodel_handle);
    2624
    … …  
    2927
    3028                if (iomodel->sub_analysis_type==HorizAnalysisEnum()){
    31                        
     29
    3230                        CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes,pmaterials, iomodel,iomodel_handle);
    3331                        CreateConstraintsDiagnosticHoriz(pconstraints,iomodel,iomodel_handle);
    3432                        CreateLoadsDiagnosticHoriz(ploads,iomodel,iomodel_handle);
    3533                        CreateParametersDiagnosticHoriz(pparameters,iomodel,iomodel_handle);
    36                                
     34
    3735                }
    3836                else if (iomodel->sub_analysis_type==VertAnalysisEnum()){
    … …  
    9189                throw ErrorException(__FUNCT__,exprintf("%s%i%s%i%s"," analysis_type: ",iomodel->analysis_type," sub_analysis_type: ",iomodel->sub_analysis_type," not supported yet!"));
    9290        }
    93                        
     91
     92        CreateParametersQmu(pparameters,iomodel,iomodel_handle);
     93
    9494}

  • issm/trunk/src/c/ModelProcessorx/Qmu/CreateParametersQmu.cpp

    r1904 r2110  
    1717void CreateParametersQmu(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){
    1818       
    19         int i;
     19        int i,j,k;
    2020       
    2121        DataSet* parameters = NULL;
    2222        Param*   param = NULL;
    2323        int      count;
     24        int      second_count;
    2425       
    2526        int*     epart=NULL;
    … …  
    3839        char* qmuerrname=NULL;
    3940        char* qmuoutname=NULL;
     41        extern int my_rank;
     42                               
     43       
     44        /*parameters for mass flux: */
     45        double* qmu_mass_flux_segments=NULL;
     46        double* my_qmu_mass_flux_segments=NULL;
     47        int num_qmu_mass_flux_segments=0;
     48        int my_num_qmu_mass_flux_segments=0;
    4049
    4150        #ifdef _SERIAL_
    … …  
    5766
    5867        if(iomodel->qmu_analysis){
     68
    5969                //name of qmu input, error and output files
    6070                qmuinname=(char*)xmalloc((strlen(iomodel->name)+strlen(".qmu.in")+1)*sizeof(char));
    … …  
    111121                #endif
    112122
     123
    113124                /*Ok, we have all the variable descriptors. Build a parameter with it: */
    114125                count++;
    … …  
    140151                #endif
    141152
     153
    142154                /*Ok, we have all the response descriptors. Build a parameter with it: */
    143155                count++;
    … …  
    156168                #endif
    157169
     170
    158171                /*partition grids in iomodel->qmu_npart parts, unless a partition is already present: */
    159172                IoModelFetchData((void**)&dpart,NULL,NULL,iomodel_handle,"part","Matrix","Mat");
    … …  
    180193                param->SetDoubleVec(dpart,iomodel->numberofnodes,1);
    181194                parameters->AddObject(param);
     195
    182196
    183197                /*Ok, now if any of the variables input from Dakota are distributed, we are going to need the parameters: */
    … …  
    204218                        }
    205219                }
     220
     221
     222                /*Deal with data needed for some responses: */
     223                for(i=0;i<iomodel->numberofresponses;i++){
     224                        char* descriptor=responsedescriptors[i];
     225                        if (strcmp(descriptor,"mass_flux")==0){
     226
     227                                /*We need the qmu_mass_flux_segments to be able to compute the mass flux: */
     228                                IoModelFetchData((void**)&qmu_mass_flux_segments,&num_qmu_mass_flux_segments,NULL,iomodel_handle,"qmu_mass_flux_segments","Matrix","Mat");
     229
     230                                #ifdef _PARALLEL_
     231                                        /*Use the element partitioning vector from the iomodel to down select qmu_mass_flux_segments to only segments that are relevant
     232                                         * to this cpu: */
     233                                        my_num_qmu_mass_flux_segments=0;
     234                                        for(j=0;j<num_qmu_mass_flux_segments;j++){
     235                                                if (  (*(qmu_mass_flux_segments+5*j+4))   == iomodel->epart[my_rank])my_qmu_mass_flux_segments++;
     236                                        }
     237                                        if(my_num_qmu_mass_flux_segments){
     238                                                my_qmu_mass_flux_segments=(double*)xcalloc(5*my_num_qmu_mass_flux_segments,sizeof(double));
     239                                                second_count=0;
     240                                                for(j=0;j<num_qmu_mass_flux_segments;j++){
     241                                                        if (*(qmu_mass_flux_segments+5*j+4)==iomodel->epart[my_rank]){
     242                                                                for(k=0;k<5;k++)*(my_qmu_mass_flux_segments+5*second_count+k)=*(qmu_mass_flux_segments+5*j+k);
     243                                                                second_count++;
     244                                                        }
     245                                                }
     246                                        }
     247
     248                                        count++;
     249                                        param= new Param(count,"qmu_mass_flux_segments",DOUBLEMAT);
     250                                        param->SetDoubleMat(my_qmu_mass_flux_segments,my_num_qmu_mass_flux_segments,5);
     251                                        parameters->AddObject(param);
     252                                #else
     253
     254                                        count++;
     255                                        param= new Param(count,"qmu_mass_flux_segments",DOUBLEMAT);
     256                                        param->SetDoubleMat(qmu_mass_flux_segments,num_qmu_mass_flux_segments,5);
     257                                        parameters->AddObject(param);
     258
     259                                #endif
     260
     261                                xfree((void**)&qmu_mass_flux_segments);
     262                                xfree((void**)&my_qmu_mass_flux_segments);
     263                        }
     264                }
     265
    206266
    207267                /*Free data: */

  • issm/trunk/src/c/io/WriteParams.cpp

    r962 r2110  
    2828        mwSize          ndim=2;
    2929        mxArray*    pfield=NULL;
     30        mxArray*    pfield2=NULL;
    3031
    3132        /*intermediary: */
    … …  
    109110                                N=param->GetN();
    110111                                pfield=mxCreateDoubleMatrix(0,0,mxREAL);
    111                                 mxSetM(pfield,M);
    112                                 mxSetN(pfield,N);
     112                                mxSetM(pfield,N);
     113                                mxSetN(pfield,M);
    113114                                mxSetPr(pfield,doublemat);
    114                                 mxSetField( dataref, 0, param->GetParameterName(),pfield);
     115                                //transpose the matrix, written directly to matlab! from C to matlab.
     116                                mexCallMATLAB(1,&pfield2, 1, &pfield, "'");
     117                                mxSetField( dataref, 0, param->GetParameterName(),pfield2);
    115118                                break;
    116119               

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

    r1881 r2110  
    5959#include "./ElementConnectivityx/ElementConnectivityx.h"
    6060#include "./OutputRiftsx/OutputRiftsx.h"
     61#include "./MassFluxx/MassFluxx.h"
    6162
    6263

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

    r1439 r2110  
    665665
    666666}
    667 
     667#undef __FUNCT__
     668#define __FUNCT__ "Beam::MassFlux"
     669double Beam::MassFlux( double* segment,double* vx,double* vy,double* vz){
     670        throw ErrorException(__FUNCT__," not supported yet!");
     671}

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

    r1188 r2110  
    8787                void  GetParameterValue(double* pvalue, double* value_list,double gauss_coord);
    8888                void  GetJacobianDeterminant(double* pJdet,double* z_list, double gauss_coord);
     89                double MassFlux(double* segment,double* vx,double* vy,double* vz);
    8990
    9091};

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

    r1188 r2110  
    4040                virtual double Misfit(void* inputs,int analysis_type,int sub_analysis_type)=0;
    4141                virtual void  ComputePressure(Vec p_g)=0;
     42                virtual double MassFlux(double* segment,double* vx,double* vy,double* vz)=0;
    4243
    4344                int           Enum();

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

    r1936 r2110  
    272272                        memcpy(marshalled_dataset,&N,sizeof(N));marshalled_dataset+=sizeof(N);
    273273                        if(M*N){
    274                                 memcpy(marshalled_dataset,doublevec,M*N*sizeof(double));
     274                                memcpy(marshalled_dataset,doublemat,M*N*sizeof(double));
    275275                                marshalled_dataset+=(M*N*sizeof(double));
    276276                        }
    … …  
    694694
    695695}
     696               
     697
     698#undef __FUNCT__
     699#define __FUNCT__ "SetDoubleMat"
     700void  Param::SetDoubleMat(double* value,int pM,int pN){
     701       
     702        if (type!=DOUBLEMAT) throw ErrorException(__FUNCT__,exprintf("%s%i"," trying to set doublematrix type",type));
     703
     704        this->M=pM;
     705        this->N=pN;
     706        if(this->M*this->N){
     707                xfree((void**)&doublemat); doublemat=(double*)xcalloc(M*N,sizeof(double));
     708                memcpy(doublemat,value,M*N*sizeof(double));
     709        }
     710}

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

    r803 r2110  
    5050                void  SetDoubleVec(double* value,int size);
    5151                void  SetDoubleVec(double* value,int size,int ndof);
     52                void  SetDoubleMat(double* value,int M,int N);
    5253                void  SetVec(Vec value);
    5354                void  SetInteger(int value);

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

    r2005 r2110  
    40284028        *phi=2*pow(epsilon_eff,2.0)*viscosity;
    40294029}
     4030
     4031
     4032#undef __FUNCT__
     4033#define __FUNCT__ "Penta::MassFlux"
     4034double Penta::MassFlux( double* segment,double* vx,double* vy,double* vz){
     4035        throw ErrorException(__FUNCT__," not supported yet!");
     4036}

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

    r1676 r2110  
    142142                void  CreatePVectorMelting( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type);
    143143                void  GetPhi(double* phi, double*  epsilon, double viscosity);
     144                double MassFlux(double* segment,double* vx,double* vy,double* vz);
    144145
    145146

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

    r1439 r2110  
    502502
    503503}
     504
     505#undef __FUNCT__
     506#define __FUNCT__ "Sing::MassFlux"
     507double Sing::MassFlux( double* segment,double* vx,double* vy,double* vz){
     508        throw ErrorException(__FUNCT__," not supported yet!");
     509}

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

    r1188 r2110  
    7979                void  GradjB(_p_Vec*, void*, int,int);
    8080                double Misfit(void*,int,int);
     81                double MassFlux(double* segment,double* vx,double* vy,double* vz);
    8182
    8283

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

    r2046 r2110  
    37093709}
    37103710
    3711 
     3711#undef __FUNCT__
     3712#define __FUNCT__ "Tria::MassFlux"
     3713double Tria::MassFlux( double* segment,double* vx,double* vy,double* vz){
     3714
     3715        int i;
     3716
     3717        const int    numgrids=3;
     3718        double mass_flux=0;
     3719        int    doflist[3];
     3720        double vx_list[3];
     3721        double vy_list[3];
     3722        double xyz_list[numgrids][3];
     3723        double gauss_1[3];
     3724        double gauss_2[3];
     3725        double normal[2];
     3726        double length;
     3727        double x1,y1,x2,y2;
     3728        double h1,h2;
     3729        double vx1,vx2,vy1,vy2;
     3730        double rho_ice;
     3731       
     3732        /*Get material parameters :*/
     3733        rho_ice=this->matpar->GetRhoIce();
     3734
     3735        /*First off, check that this segment belongs to this element: */
     3736        if ((int)*(segment+4)!=this->id)throw ErrorException(__FUNCT__,exprintf("%s%i%s%i","error message: segment with id ",(int)*(segment+4)," does not belong to element with id:",this->id));
     3737
     3738        /*Recover segment node locations: */
     3739        x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
     3740       
     3741        /*Get xyz list: */
     3742        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     3743
     3744        /*recover velocity at three element nodes: */
     3745        this->GetDofList1(&doflist[0]);
     3746        for(i=0;i<3;i++){
     3747                vx_list[i]=vx[doflist[i]];
     3748                vy_list[i]=vy[doflist[i]];
     3749        }
     3750
     3751        /*get area coordinates of 0 and 1 locations: */
     3752        for(i=0;i<3;i++){
     3753                gauss_1[i]=this->GetAreaCoordinate(x1,y1,i+1);
     3754                gauss_2[i]=this->GetAreaCoordinate(x2,y2,i+1);
     3755        }
     3756
     3757        /*get normal of segment: */
     3758        normal[0]=cos(atan2(x1-x2,y2-y1));
     3759        normal[1]=sin(atan2(x1-x2,y2-y1));
     3760
     3761        /*get length of segment: */
     3762        length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2));
     3763
     3764        /*get thickness and velocity at two segment extremities: */
     3765        GetParameterValue(&h1, &h[0],gauss_1);
     3766        GetParameterValue(&h2, &h[0],gauss_2);
     3767        GetParameterValue(&vx1, &vx_list[0],gauss_1);
     3768        GetParameterValue(&vy1, &vy_list[0],gauss_1);
     3769        GetParameterValue(&vx2, &vx_list[0],gauss_2);
     3770        GetParameterValue(&vy2, &vy_list[0],gauss_2);
     3771
     3772
     3773        mass_flux= rho_ice*length*( 
     3774                                  (1.0/3.0*(h1-h2)*(vx1-vx2)+1.0/2.0*h2*(vx1-vx2)+1.0/2.0*(h1-h2)*vx2+h2*vx2)*normal[0]+
     3775                                  (1.0/3.0*(h1-h2)*(vy1-vy2)+1.0/2.0*h2*(vy1-vy2)+1.0/2.0*(h1-h2)*vy2+h2*vy2)*normal[1]
     3776                                );
     3777
     3778        return mass_flux;
     3779}
     3780
     3781#undef __FUNCT__
     3782#define __FUNCT__ "Tria::GetArea"
     3783double Tria::GetArea(void){
     3784
     3785        double area=0;
     3786        const int    numgrids=3;
     3787        double xyz_list[numgrids][3];
     3788        double x1,y1,x2,y2,x3,y3;
     3789
     3790        /*Get xyz list: */
     3791        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     3792        x1=xyz_list[0][0]; y1=xyz_list[0][1];
     3793        x2=xyz_list[1][0]; y2=xyz_list[1][1];
     3794        x3=xyz_list[2][0]; y3=xyz_list[2][1];
     3795 
     3796        return x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1;
     3797}
     3798
     3799#undef __FUNCT__
     3800#define __FUNCT__ "Tria::GetAreaCoordinate"
     3801double Tria::GetAreaCoordinate(double x, double y, int which_one){
     3802
     3803        double area=0;
     3804        const int    numgrids=3;
     3805        double xyz_list[numgrids][3];
     3806        double x1,y1,x2,y2,x3,y3;
     3807
     3808        /*Get area: */
     3809        area=this->GetArea();
     3810
     3811        /*Get xyz list: */
     3812        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     3813        x1=xyz_list[0][0]; y1=xyz_list[0][1];
     3814        x2=xyz_list[1][0]; y2=xyz_list[1][1];
     3815        x3=xyz_list[2][0]; y3=xyz_list[2][1];
     3816
     3817        if(which_one==1){
     3818                /*Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area*/
     3819                return ((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area;
     3820        }
     3821        else if(which_one==2){
     3822                /*Get second area coordinate = det(x1-x3  x-x3 ; y1-y3   y-y3)/area*/
     3823                return ((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area;
     3824        }
     3825        else if(which_one==3){
     3826                /*Get third  area coordinate 1-area1-area2: */
     3827                return 1-((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area -((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area;
     3828        }
     3829        else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n"," error message: area coordinate ",which_one," done not exist!"));
     3830}
     3831
     3832

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

    r1211 r2110  
    122122                void  CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    123123                void  CreatePVectorPrognostic(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type);
     124                double MassFlux(double* segment,double* vx,double* vy,double* vz);
     125                double GetArea(void);
     126                double GetAreaCoordinate(double x, double y, int which_one);
    124127
    125128

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

    r2016 r2110  
    3636        Vec slopey=NULL;
    3737        Vec riftproperties=NULL;
     38        double* u_g_initial=NULL;
    3839
    3940        /*flags: */
    4041        int debug=0;
     42        int qmu_analysis=0;
    4143        int dim=-1;
    4244        int ishutter=0;
    … …  
    6870        model->FindParam(&stokesreconditioning,"stokesreconditioning");
    6971        model->FindParam(&numrifts,"numrifts");
     72        model->FindParam(&qmu_analysis,"qmu_analysis");
    7073
    7174        /*recover fem models: */
    … …  
    8083        fem_sl->FindParam((void*)&numberofdofspernode_sl,"numberofdofspernode");
    8184        fem_ds->FindParam((void*)&numberofdofspernode_ds,"numberofdofspernode");
     85
     86        //for qmu analysis, be sure the velocity input we are starting from  is the one in the parameters: */
     87        if(qmu_analysis){
     88                model->FindParam(&u_g_initial,"u_g",DiagnosticAnalysisEnum(),HorizAnalysisEnum());
     89                inputs->Add("velocity",u_g_initial,3,numberofnodes);
     90        }
    8291
    8392        if(ishutter){
    … …  
    202211        VecFree(&pg);
    203212        xfree((void**)&dofset);
     213        xfree((void**)&u_g_initial);
    204214
    205215}

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

    r2021 r2110  
    274274        md.variabledescriptors=NaN;
    275275        md.responsedescriptors=NaN;
     276        md.qmu_mass_flux_profile=NaN;
     277        md.qmu_mass_flux_segments=NaN;
    276278
    277279        %Ice solver string

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

    r213 r2110  
    3333        fwrite(fid,s(2),'int');
    3434        if s(1)*s(2),
    35                 fwrite(fid,data','double'); %get to the "c" convenction, hence the transpose
     35                fwrite(fid,data','double'); %get to the "c" convention, hence the transpose
    3636        end
    3737elseif strcmpi(data_type,'Integer'),

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

    r2049 r2110  
    5656        factors=importancefactors(npart)';
    5757end
     58
     59if numel(factors)==md.numberofgrids,
     60        %get areas for each vertex.
     61        aire=GetAreas(md.elements,md.x,md.y);
     62        num_elements_by_node=md.nodeconnectivity(:,end);
     63        grid_aire=zeros(md.numberofgrids,1);
     64        for i=1:md.numberofgrids,
     65                for j=1:num_elements_by_node(i),
     66                        grid_aire(i)=grid_aire(i)+aire(md.nodeconnectivity(i,j));
     67                end
     68        end
     69        factors=factors./grid_aire;
     70end

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

    r1805 r2110  
    2929        count=count+numpairsforthisrift;
    3030end
     31
     32

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

    r832 r2110  
    77
    88function  m=CreateFEMModel(md)
    9 
     9       
    1010        displaystring(md.debug,'\n   reading data from model %s...',md.name);
    1111        [m.elements,m.nodes,m.constraints,m.loads,m.materials,parameters]=ModelProcessor(md);
    12        
     12
    1313        displaystring(md.debug,'%s','   generating degrees of freedom...');
    1414        [m.nodes,m.part,m.tpart]=Dof(m.elements,m.nodes,parameters);

  • issm/trunk/src/m/solutions/dakota/preqmu.m

    r1029 r2110  
    114114md.responsedescriptors=responsedescriptors;
    115115
     116%now, we have to provide all the info necessary for the solutions to compute the responses. For ex, if mass_flux
     117%is a response, we need a profile of points. For max_vel, we don't need anything.
     118md=process_qmu_response_data(md);

  • issm/trunk/src/m/solutions/dakota/qmumarshall.m

    r1109 r2110  
    6969end
    7070
     71%write response specific data
     72count=0;
     73for i=1:length(response_fieldnames),
     74        field_name=response_fieldnames{i};
     75        fieldresponses=responses.(field_name);
     76        for j=1:length(fieldresponses),
     77                descriptor=fieldresponses(j).descriptor;
     78                if strcmpi(descriptor,'mass_flux'),
     79                        WriteData(fid,md.qmu_mass_flux_segments,'Mat','qmu_mass_flux_segments');
     80                end
     81        end
     82end
     83
    7184%write part and npart to disk
    7285WriteData(fid,md.npart,'Integer','npart');

  • issm/trunk/src/m/solutions/dakota/qmuresponse.m

    r963 r2110  
    44if strcmpi(descriptor,'max_vel'),
    55        response=max(results.vel);
     6elseif strcmpi(descriptor,'min_vel'),
     7        response=min(results.vel);
     8elseif strcmpi(descriptor,'max_vx'),
     9        response=max(results.vx);
     10elseif strcmpi(descriptor,'min_vx'),
     11        response=min(results.vx);
     12elseif strcmpi(descriptor,'max_vy'),
     13        response=max(results.vy);
     14elseif strcmpi(descriptor,'min_vy'),
     15        response=min(results.vy);
     16elseif strcmpi(descriptor,'mass_flux'),
     17        %call mass flux module.
     18        m_dh=models.dh;
     19        m_dhu=models.dhu;
     20        m_ds=models.ds;
     21        ishutter=m_dh.parameters.ishutter;
     22        ismacayealpattyn=m_dhu.parameters.ismacayealpattyn;
     23        isstokes=m_ds.parameters.isstokes;
     24        if ishutter,
     25                response=MassFlux(m_dhu.elements,m_dhu.nodes,m_dhu.loads,m_dhu.materials,m_dhu.parameters,results);
     26        elseif ismacayealpattyn,
     27                response=MassFlux(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,m_dh.parameters,results);
     28        elseif isstokes,
     29                response=MassFlux(m_ds.elements,m_ds.nodes,m_ds.loads,m_ds.materials,m_ds.parameters,results);
     30        else
     31                error('qmuresponse error message: unsupported analysis type for mass_flux computation!');
     32        end
    633else
    734        error(['qmuresponse error message: unknow descriptor ' descriptor]);

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

    r2047 r2110  
    2222                                InterpFromMesh2d \
    2323                                InterpFromMesh3d \
     24                                MassFlux\
    2425                                Mergesolutionfromftog\
    2526                                MeshPartition\
    … …  
    137138                          Mergesolutionfromftog/Mergesolutionfromftog.h
    138139
     140MassFlux_SOURCES = MassFlux/MassFlux.cpp\
     141                          MassFlux/MassFlux.h
     142
    139143MeshPartition_SOURCES = MeshPartition/MeshPartition.cpp\
    140144                          MeshPartition/MeshPartition.h

  • issm/trunk/test/Verification/IceSheetIceFrontM2dDakota_25/Square.par

    r2056 r2110  
    6060md.qmu_analysis=1;
    6161
    62 md.eps_rel=10^-15; %tighten for qmu analysese
     62md.eps_rel=10^-10; %tighten for qmu analysese

  • issm/trunk/test/Verification/IceSheetIceFrontM2dDakota_25/runme.m

    r2056 r2110  
    6060                %initialize model
    6161                md=model;
    62                 md=mesh(md,'DomainOutline.exp',100000);
     62                md=mesh(md,'DomainOutline.exp',60000);
    6363                md=geography(md,'','');
    6464                md=parameterize(md,'Square.par');
    … …  
    7575                md=tres(md,'dakota');
    7676                md.results.dakota.importancefactors=importancefactors(md,'drag','max_vel');
     77                varargout{1}=md;
     78                return;
    7779
    7880                %compute fields to be checked
Note: See TracChangeset for help on using the changeset viewer.