Ignore:
Timestamp:
11/04/16 13:48:43 (8 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 21337

Location:
issm/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c

    • Property svn:ignore

      •  

        old new  
        55*.obj
        66*.exe
        7 issm
        8 kriging
        97appscan.*
        108ouncemake*
        … …  
        1917*.dirstamp
        2018.libs
         19issm
         20kriging
         21issm_slr

  • issm/trunk/src/c/classes/Elements/Penta.cpp

    r20500 r21341  
    3030
    3131        /*Checks in debugging mode*/
    32         _assert_(iomodel->Data(MeshUpperelementsEnum));
    33         _assert_(iomodel->Data(MeshLowerelementsEnum));
     32        _assert_(iomodel->Data("md.mesh.upperelements"));
     33        _assert_(iomodel->Data("md.mesh.lowerelements"));
    3434
    3535        /*id: */
    … …  
    3838
    3939        /*Build neighbors list*/
    40         if (xIsNan<IssmDouble>(iomodel->Data(MeshUpperelementsEnum)[index]) || iomodel->Data(MeshUpperelementsEnum)[index]==-1.) penta_elements_ids[1]=this->id; //upper penta is the same penta
    41         else                                    penta_elements_ids[1]=reCast<int,IssmDouble>((iomodel->Data(MeshUpperelementsEnum)[index]));
    42         if (xIsNan<IssmDouble>(iomodel->Data(MeshLowerelementsEnum)[index]) || iomodel->Data(MeshLowerelementsEnum)[index]==-1.) penta_elements_ids[0]=this->id; //lower penta is the same penta
    43         else                                    penta_elements_ids[0]=reCast<int,IssmDouble>((iomodel->Data(MeshLowerelementsEnum)[index]));
     40        if (xIsNan<IssmDouble>(iomodel->Data("md.mesh.upperelements")[index]) || iomodel->Data("md.mesh.upperelements")[index]==-1.) penta_elements_ids[1]=this->id; //upper penta is the same penta
     41        else                                    penta_elements_ids[1]=reCast<int,IssmDouble>((iomodel->Data("md.mesh.upperelements")[index]));
     42        if (xIsNan<IssmDouble>(iomodel->Data("md.mesh.lowerelements")[index]) || iomodel->Data("md.mesh.lowerelements")[index]==-1.) penta_elements_ids[0]=this->id; //lower penta is the same penta
     43        else                                    penta_elements_ids[0]=reCast<int,IssmDouble>((iomodel->Data("md.mesh.lowerelements")[index]));
    4444        this->InitHookNeighbors(penta_elements_ids);
    4545
    … …  
    289289                /*Compute strain rate viscosity and pressure: */
    290290                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    291                 this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     291                this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    292292                pressure_input->GetInputValue(&pressure,gauss);
    293293
    … …  
    329329        IssmDouble      tau_xz[NUMVERTICES];
    330330        IssmDouble      tau_yz[NUMVERTICES];
     331        IssmDouble      tau_eff[NUMVERTICES];
    331332        GaussPenta* gauss=NULL;
    332333
    … …  
    346347                /*Compute strain rate viscosity and pressure: */
    347348                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    348                 this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     349                this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    349350
    350351                /*Compute Stress*/
    … …  
    355356                tau_xz[iv]=2*viscosity*epsilon[4];
    356357                tau_yz[iv]=2*viscosity*epsilon[5];
     358
     359                tau_eff[iv] = tau_xx[iv]*tau_xx[iv] + tau_yy[iv]*tau_yy[iv] + tau_zz[iv]*tau_zz[iv] +
     360                  2*tau_xy[iv]*tau_xy[iv] + 2*tau_xz[iv]*tau_xz[iv] + 2*tau_yz[iv]*tau_yz[iv];
     361
     362                tau_eff[iv] = sqrt(tau_eff[iv]/2.);
    357363        }
    358364
    … …  
    364370        this->inputs->AddInput(new PentaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum));
    365371        this->inputs->AddInput(new PentaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
     372        this->inputs->AddInput(new PentaInput(DeviatoricStresseffectiveEnum,&tau_eff[0],P1Enum));
    366373
    367374        /*Clean up and return*/
    … …  
    398405                /*Compute strain rate viscosity and pressure: */
    399406                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    400                 this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     407                this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    401408                pressure_input->GetInputValue(&pressure,gauss);
    402409
    … …  
    640647                        /*Check for basal force only if grounded and touching GL*/
    641648                        this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
    642                         this->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
     649                        this->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
    643650                        pressure_input->GetInputValue(&pressure, gauss);
    644651                        base_input->GetInputValue(&base, gauss); _assert_(base<0.);
    … …  
    10731080
    10741081        /*Get out if this is not an element input*/
    1075         if(!IsInput(control_enum)) return;
     1082        if(!IsInput(control_enum)) _error_("Enum "<<EnumToStringx(control_enum)<<" is not in IsInput");
    10761083
    10771084        /*Prepare index list*/
    … …  
    13961403        IssmDouble  yts;
    13971404        bool    control_analysis;
     1405        char**  controls = NULL;
    13981406        int     num_control_type,num_responses;
    13991407
    14001408        /*Fetch parameters: */
    1401         iomodel->Constant(&yts,ConstantsYtsEnum);
    1402         iomodel->Constant(&control_analysis,InversionIscontrolEnum);
    1403         if(control_analysis) iomodel->Constant(&num_control_type,InversionNumControlParametersEnum);
    1404         if(control_analysis) iomodel->Constant(&num_responses,InversionNumCostFunctionsEnum);
     1409        iomodel->FindConstant(&yts,"md.constants.yts");
     1410        iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
     1411        if(control_analysis) iomodel->FindConstant(&num_control_type,"md.inversion.num_control_parameters");
     1412        if(control_analysis) iomodel->FindConstant(&num_responses,"md.inversion.num_cost_functions");
    14051413
    14061414        /*Recover vertices ids needed to initialize inputs*/
    … …  
    14111419
    14121420        /*Control Inputs*/
    1413         if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
     1421        if (control_analysis){
     1422                iomodel->FindConstant(&controls,NULL,"md.inversion.control_parameters");
    14141423                for(i=0;i<num_control_type;i++){
    1415                         int control = reCast<int>(iomodel->Data(InversionControlParametersEnum)[i]);
     1424                        _assert_(controls[i]);
     1425                        int control = StringToEnumx(controls[i]);
    14161426                        switch(control){
    1417                                 /*yts conversion*/
    14181427                                case BalancethicknessThickeningRateEnum:
    1419                                 case VxEnum:
    1420                                 case VyEnum:
    1421                                         if(iomodel->Data(control)){
    1422                                                 for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(control)[penta_vertex_ids[j]-1];
    1423                                                 for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
    1424                                                 for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
    1425                                                 this->inputs->AddInput(new ControlInput(control,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1428                                        if (iomodel->Data("md.balancethickness.thickening_rate")){
     1429                                                for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data("md.balancethickness.thickening_rate")[penta_vertex_ids[j]-1];
     1430                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
     1431                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
     1432                                                this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    14261433                                        }
    14271434                                        break;
    1428 
    1429                                 /*No yts conversion*/
    1430                                 case ThicknessEnum:
    1431                                 case FrictionCoefficientEnum:
    1432                           case FrictionAsEnum: 
    1433                                 case MaterialsRheologyBEnum:
    1434                                         if(iomodel->Data(control)){
    1435                                                 for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(control)[penta_vertex_ids[j]-1];
    1436                                                 for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
    1437                                                 for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
    1438                                                 this->inputs->AddInput(new ControlInput(control,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1435                                case VxEnum:
     1436                                        if (iomodel->Data("md.initialization.vx")){
     1437                                                for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data("md.initialization.vx")[penta_vertex_ids[j]-1];
     1438                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
     1439                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
     1440                                                this->inputs->AddInput(new ControlInput(VxEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    14391441                                        }
    14401442                                        break;
    1441 
    1442                                 /*Special cases (depth averaged quantities)*/
     1443                                case VyEnum:
     1444                                        if (iomodel->Data("md.initialization.vy")){
     1445                                                for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data("md.initialization.vy")[penta_vertex_ids[j]-1];
     1446                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
     1447                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i]/yts;
     1448                                                this->inputs->AddInput(new ControlInput(VyEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1449                                        }
     1450                                        break;
     1451                                case FrictionCoefficientEnum:
     1452                                        if (iomodel->Data("md.friction.coefficient")){
     1453                                                for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data("md.friction.coefficient")[penta_vertex_ids[j]-1];
     1454                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i];
     1455                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i];
     1456                                                this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
     1457                                        }
     1458                                        break;
    14431459                                case MaterialsRheologyBbarEnum:
    1444                                         if(iomodel->Data(MaterialsRheologyBEnum)){
    1445                                                 for(j=0;j<NUMVERTICES;j++) nodeinputs[j]=iomodel->Data(MaterialsRheologyBEnum)[penta_vertex_ids[j]-1];
    1446                                                 for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
    1447                                                 for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
     1460                                        if(iomodel->Data("md.materials.rheology_B")){
     1461                                                for(j=0;j<NUMVERTICES;j++) nodeinputs[j]=iomodel->Data("md.materials.rheology_B")[penta_vertex_ids[j]-1];
     1462                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i];
     1463                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i];
    14481464                                                this->inputs->AddInput(new ControlInput(MaterialsRheologyBEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    14491465                                        }
    14501466                                        break;
    14511467                                case DamageDbarEnum:
    1452                                         if(iomodel->Data(DamageDEnum)){
    1453                                                 for(j=0;j<NUMVERTICES;j++) nodeinputs[j]=iomodel->Data(DamageDEnum)[penta_vertex_ids[j]-1];
    1454                                                 for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
    1455                                                 for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(penta_vertex_ids[j]-1)*num_control_type+i];
     1468                                        if(iomodel->Data("md.damage.D")){
     1469                                                for(j=0;j<NUMVERTICES;j++) nodeinputs[j]=iomodel->Data("md.damage.D")[penta_vertex_ids[j]-1];
     1470                                                for(j=0;j<NUMVERTICES;j++)cmmininputs[j]=iomodel->Data("md.inversion.min_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i];
     1471                                                for(j=0;j<NUMVERTICES;j++)cmmaxinputs[j]=iomodel->Data("md.inversion.max_parameters")[(penta_vertex_ids[j]-1)*num_control_type+i];
    14561472                                                this->inputs->AddInput(new ControlInput(DamageDEnum,PentaInputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
    14571473                                        }
    … …  
    14611477                        }
    14621478                }
     1479                for(i=0;i<num_control_type;i++) xDelete<char>(controls[i]);
     1480                xDelete<char*>(controls);
    14631481        }
    14641482
    14651483        /*Need to know the type of approximation for this element*/
    1466         if(iomodel->Data(FlowequationElementEquationEnum)){
    1467                 this->inputs->AddInput(new IntInput(ApproximationEnum,reCast<int>(iomodel->Data(FlowequationElementEquationEnum)[index])));
     1484        if(iomodel->Data("md.flowequation.element_equation")){
     1485                this->inputs->AddInput(new IntInput(ApproximationEnum,IoCodeToEnumElementEquation(reCast<int>(iomodel->Data("md.flowequation.element_equation")[index]))));
    14681486        }
    14691487
    14701488        /*DatasetInputs*/
    1471         if (control_analysis && iomodel->Data(InversionCostFunctionsCoefficientsEnum)) {
     1489        if(control_analysis && iomodel->Data("md.inversion.cost_functions_coefficients")) {
     1490
     1491                /*Generate cost functions associated with the iomodel*/
     1492                char**  cost_functions                  = NULL;
     1493                int*            cost_functions_enums = NULL;
     1494                int             num_cost_functions;
     1495
     1496                iomodel->FindConstant(&num_cost_functions,"md.inversion.num_cost_functions");
     1497                iomodel->FindConstant(&cost_functions,&num_cost_functions,"md.inversion.cost_functions");
     1498                if(num_cost_functions<1) _error_("No cost functions found");
     1499                cost_functions_enums=xNew<int>(num_cost_functions);
     1500                for(j=0;j<num_cost_functions;j++){ cost_functions_enums[j]=StringToEnumx(cost_functions[j]); }
    14721501
    14731502                /*Create inputs and add to DataSetInput*/
    14741503                DatasetInput* datasetinput=new DatasetInput(InversionCostFunctionsCoefficientsEnum);
    14751504                for(i=0;i<num_responses;i++){
    1476                         for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data(InversionCostFunctionsCoefficientsEnum)[(penta_vertex_ids[j]-1)*num_responses+i];
    1477                         datasetinput->AddInput(new PentaInput(InversionCostFunctionsCoefficientsEnum,nodeinputs,P1Enum),reCast<int>(iomodel->Data(InversionCostFunctionsEnum)[i]));
     1505                        for(j=0;j<NUMVERTICES;j++)nodeinputs[j]=iomodel->Data("md.inversion.cost_functions_coefficients")[(penta_vertex_ids[j]-1)*num_responses+i];
     1506                        datasetinput->AddInput(new PentaInput(InversionCostFunctionsCoefficientsEnum,nodeinputs,P1Enum),cost_functions_enums[i]);
    14781507                }
    14791508
    14801509                /*Add datasetinput to element inputs*/
    14811510                this->inputs->AddInput(datasetinput);
     1511
     1512                /*Free resources*/
     1513                for(int j=0;j<num_cost_functions;j++) xDelete<char>(cost_functions[j]);
     1514                xDelete<char*>(cost_functions);
     1515                xDelete<int>(cost_functions_enums);
    14821516        }
    14831517}
    … …  
    14991533                values[i]=solution[doflist[i]];
    15001534                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
     1535                if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
    15011536        }
    15021537
    … …  
    15291564                values[i+numdof2d]=values[i];
    15301565                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
     1566                if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
    15311567        }
    15321568
    … …  
    15551591
    15561592        /*Check that name is an element input*/
    1557         if (!IsInput(name)) return;
     1593        if(!IsInput(name)) _error_("Enum "<<EnumToStringx(name)<<" is not in IsInput");
    15581594
    15591595        switch(type){
    … …  
    15821618                                values[i]=vector[doflist[i]];
    15831619                                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
     1620                                if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
    15841621                        }
    15851622                        /*Add input to the element: */
    … …  
    15941631                                values[i]=vector[nodes[i]->Sid()];
    15951632                                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
     1633                                if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
    15961634                        }
    15971635                        /*Add input to the element: */
    … …  
    22452283        _assert_(this->IsOnBase());
    22462284
    2247         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
    2248         if(this->material->IsDamage())this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
     2285        switch(this->material->ObjectEnum()){
     2286                case MaticeEnum:
     2287                        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
     2288                        if(this->material->IsDamage())this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
     2289                        break;
     2290                case MatestarEnum:
     2291                        this->InputDepthAverageAtBase(MaterialsRheologyKoEnum,MaterialsRheologyKobarEnum);
     2292                        this->InputDepthAverageAtBase(MaterialsRheologyEcEnum,MaterialsRheologyEcbarEnum);
     2293                        this->InputDepthAverageAtBase(MaterialsRheologyEsEnum,MaterialsRheologyEsbarEnum);
     2294                        break;
     2295                default:
     2296                        _error_("not supported yet");
     2297        }
    22492298        if(this->inputs->GetInput(VxEnum)) this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
    22502299        if(this->inputs->GetInput(VyEnum)) this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
    … …  
    22522301        if(this->inputs->GetInput(CalvingrateyEnum)) this->InputDepthAverageAtBase(CalvingrateyEnum,CalvingrateyAverageEnum);
    22532302        Tria* tria=(Tria*)SpawnTria(0,1,2);
    2254         this->inputs->DeleteInput(MaterialsRheologyBbarEnum);
    2255         this->inputs->DeleteInput(DamageDbarEnum);
     2303        switch(this->material->ObjectEnum()){
     2304                case MaticeEnum:
     2305                        this->inputs->DeleteInput(MaterialsRheologyBbarEnum);
     2306                        this->inputs->DeleteInput(DamageDbarEnum);
     2307                        break;
     2308                case MatestarEnum:
     2309                        break;
     2310                default:
     2311                        _error_("not supported yet");
     2312        }
    22562313        this->inputs->DeleteInput(VxAverageEnum);
    22572314        this->inputs->DeleteInput(VyAverageEnum);
    … …  
    26802737
    26812738        /*Fetch parameters: */
    2682         iomodel->Constant(&yts,ConstantsYtsEnum);
    2683         iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
     2739        iomodel->FindConstant(&yts,"md.constants.yts");
     2740        iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota");
    26842741
    26852742        /*Checks if debuging*/
    … …  
    30343091
    30353092                case StressbalanceAnalysisEnum:
    3036                         _assert_(iomodel->Data(FlowequationElementEquationEnum));
    3037 
    3038                         if(*(iomodel->Data(FlowequationElementEquationEnum)+index)==HOFSApproximationEnum){
     3093                        _assert_(iomodel->Data("md.flowequation.element_equation"));
     3094
     3095                        if((IoCodeToEnumElementEquation(reCast<int>(iomodel->Data("md.flowequation.element_equation")[index])))==HOFSApproximationEnum){
    30393096                                /*Create VzHO and VzFS Enums*/
    3040                                 if(iomodel->Data(VzEnum) && iomodel->Data(FlowequationBorderFSEnum)){
    3041                                         for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*iomodel->Data(FlowequationBorderFSEnum)[penta_vertex_ids[i]-1];
     3097                                if(iomodel->Data("md.initialization.vz") && iomodel->Data("md.flowequation.borderFS")){
     3098                                        for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data("md.initialization.vz")[penta_vertex_ids[i]-1]*iomodel->Data("md.flowequation.borderFS")[penta_vertex_ids[i]-1];
    30423099                                        this->inputs->AddInput(new PentaInput(VzFSEnum,nodeinputs,P1Enum));
    3043                                         for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*(1-iomodel->Data(FlowequationBorderFSEnum)[penta_vertex_ids[i]-1]);
     3100                                        for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data("md.initialization.vz")[penta_vertex_ids[i]-1]*(1-iomodel->Data("md.flowequation.borderFS")[penta_vertex_ids[i]-1]);
    30443101                                        this->inputs->AddInput(new PentaInput(VzHOEnum,nodeinputs,P1Enum));
    30453102                                }
    … …  
    30503107                                }
    30513108                        }
    3052                         if(*(iomodel->Data(FlowequationElementEquationEnum)+index)==SSAFSApproximationEnum){
     3109                        if((IoCodeToEnumElementEquation(reCast<int>(iomodel->Data("md.flowequation.element_equation")[index])))==SSAFSApproximationEnum){
    30533110                                /*Create VzSSA and VzFS Enums*/
    3054                                 if(iomodel->Data(VzEnum) && iomodel->Data(FlowequationBorderFSEnum)){
    3055                                         for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*iomodel->Data(FlowequationBorderFSEnum)[penta_vertex_ids[i]-1];
     3111                                if(iomodel->Data("md.initialization.vz") && iomodel->Data("md.flowequation.borderFS")){
     3112                                        for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data("md.initialization.vz")[penta_vertex_ids[i]-1]*iomodel->Data("md.flowequation.borderFS")[penta_vertex_ids[i]-1];
    30563113                                        this->inputs->AddInput(new PentaInput(VzFSEnum,nodeinputs,P1Enum));
    3057                                         for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*(1-iomodel->Data(FlowequationBorderFSEnum)[penta_vertex_ids[i]-1]);
     3114                                        for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data("md.initialization.vz")[penta_vertex_ids[i]-1]*(1-iomodel->Data("md.flowequation.borderFS")[penta_vertex_ids[i]-1]);
    30583115                                        this->inputs->AddInput(new PentaInput(VzSSAEnum,nodeinputs,P1Enum));
    30593116                                }
    … …  
    31743231        _assert_(gauss->Enum()==GaussPentaEnum);
    31753232        this->StrainRateFS(&epsilon[0],xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input);
    3176         this->ViscosityFS(&viscosity,3,xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input);
     3233        this->material->ViscosityFS(&viscosity,3,xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input);
    31773234        GetPhi(&phi,&epsilon[0],viscosity);
    31783235
    … …  
    33453402
    33463403#ifdef _HAVE_DAKOTA_
     3404void       Penta::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
     3405
     3406        int             i,t,row;
     3407        IssmDouble      time;
     3408        TransientInput *transientinput = NULL;
     3409        IssmDouble      values[6];
     3410
     3411        /*Check that name is an element input*/
     3412        if(!IsInput(name)) _error_("Enum "<<EnumToStringx(name)<<" is not in IsInput");
     3413
     3414        switch(type){
     3415
     3416                case VertexEnum:
     3417                        /*Create transient input: */
     3418                        for(t=0;t<ncols;t++){ //ncols is the number of times
     3419
     3420                                /*create input values: */
     3421                                for(i=0;i<6;i++){
     3422                                        row=this->vertices[i]->Sid();
     3423                                        values[i]=matrix[ncols*row+t];
     3424                                }
     3425
     3426                                /*time:*/
     3427                                time=matrix[(nrows-1)*ncols+t];
     3428
     3429                                if(t==0) transientinput=new TransientInput(name);
     3430                                transientinput->AddTimeInput(new PentaInput(name,values,P1Enum),time);
     3431                                transientinput->Configure(parameters);
     3432                        }
     3433                        this->inputs->AddInput(transientinput);
     3434                        break;
     3435
     3436                default:
     3437                        _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
     3438        }
     3439
     3440}
     3441/*}}}*/
    33473442void       Penta::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/
    33483443
    … …  
    33503445
    33513446        /*Check that name is an element input*/
    3352         if (!IsInput(name)) return;
     3447        if(!IsInput(name)) _error_("Enum "<<EnumToStringx(name)<<" is not in IsInput");
    33533448
    33543449        switch(type){
    … …  
    34453540}
    34463541/*}}}*/
    3447 void       Penta::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
    3448 
    3449         int             i,t,row;
    3450         IssmDouble      time;
    3451         TransientInput *transientinput = NULL;
    3452         IssmDouble      values[6];
    3453 
    3454         /*Check that name is an element input*/
    3455         if (!IsInput(name)) return;
    3456 
    3457         switch(type){
    3458 
    3459                 case VertexEnum:
    3460                         /*Create transient input: */
    3461                         for(t=0;t<ncols;t++){ //ncols is the number of times
    3462 
    3463                                 /*create input values: */
    3464                                 for(i=0;i<6;i++){
    3465                                         row=this->vertices[i]->Sid();
    3466                                         values[i]=matrix[ncols*row+t];
    3467                                 }
    3468 
    3469                                 /*time:*/
    3470                                 time=matrix[(nrows-1)*ncols+t];
    3471 
    3472                                 if(t==0) transientinput=new TransientInput(name);
    3473                                 transientinput->AddTimeInput(new PentaInput(name,values,P1Enum),time);
    3474                                 transientinput->Configure(parameters);
    3475                         }
    3476                         this->inputs->AddInput(transientinput);
    3477                         break;
    3478 
    3479                 default:
    3480                         _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
    3481         }
    3482 
    3483 }
    3484 /*}}}*/
    34853542#endif
Note: See TracChangeset for help on using the changeset viewer.