source: issm/oecreview/Archive/21724-22754/ISSM-21930-21931.diff

Last change on this file was 22755, checked in by Mathieu Morlighem, 7 years ago

CHG: added 21724-22754
File size: 32.2 KB

  • ../trunk-jpl/src/c/bamg/GeomEdge.h

     
    2626
    2727                        //Methods
    2828                        R2     F(double theta) const ; // parametrization of the curve edge
    29                         double R1tg(double theta,R2 &t) const ; // 1/radius of curvature + tangente
    3029                        int    Cracked() const;
    3130                        int    TgA()     const;
    3231                        int    TgB()     const;

  • ../trunk-jpl/src/c/bamg/GeomEdge.cpp

     
    6262        int    GeomEdge::Mark()    const  {/*{{{*/
    6363                return type &16;
    6464        }/*}}}*/
    65         double GeomEdge::R1tg(double theta,R2 & t) const{/*{{{*/
    66                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/R1tg)*/
    67                 // 1/R of radius of cuvature
    68 
    69                 R2 A=v[0]->r,B=v[1]->r;
    70                 double dca,dcb,dcta,dctb;
    71                 double ddca,ddcb,ddcta,ddctb;
    72                 double tt = theta*theta;
    73 
    74                 //check theta
    75                 _assert_(theta>=0 && theta<=1);
    76 
    77                 if (TgA()){
    78                         if (TgB()){
    79                                 // Tangent A and B provided:
    80                                 // interpolation d'hermite
    81                                 dcb = 6*theta*(1-theta);
    82                                 ddcb = 6*(1-2*theta);
    83                                 dca = -dcb;
    84                                 ddca = -ddcb;
    85                                 dcta =  (3*theta - 4)*theta + 1;
    86                                 ddcta=6*theta-4;
    87                                 dctb = 3*tt - 2*theta;
    88                                 ddctb = 6*theta-2;
    89                         }
    90                         else {
    91                                 //Tangent A provided but tangent B not provided
    92                                 // 1-t*t, t-t*t, t*t
    93                                 double t = theta;
    94                                 dcb = 2*t;
    95                                 ddcb = 2;
    96                                 dca = -dcb;
    97                                 ddca = -2;
    98                                 dcta = 1-dcb;
    99                                 ddcta = -ddcb;
    100                                 dctb=0;   
    101                                 ddctb=0;   
    102                         }
    103                 }
    104                 else{
    105                         if (TgB()){
    106                                 //Tangent B provided but tangent A not provided
    107                                 double t = 1-theta;
    108                                 dca = -2*t;
    109                                 ddca = 2;
    110                                 dcb = -dca;
    111                                 ddcb = -2;
    112                                 dctb = 1+dca;
    113                                 ddctb= ddca;
    114                                 dcta =0;
    115                                 ddcta =0;
    116                         }
    117                         else {
    118                                 //Neither thangent A nor tangent B provided
    119                                 // lagrange P1
    120                                 t=B-A;
    121                                 return 0;
    122                         }
    123                 }
    124                 R2 d  =  A*dca  + B*dcb  + tg[0]* dcta  + tg[1] * dctb;
    125                 R2 dd =  A*ddca + B*ddcb + tg[0]* ddcta + tg[1] * ddctb;
    126                 double d2=(d,d);
    127                 double sd2 = sqrt(d2);
    128                 t=d;
    129                 if(d2>1.0e-20){
    130                         t/=sd2;
    131                         return Abs(Det(d,dd))/(d2*sd2);
    132                 }
    133                 else return 0;
    134         }
    135         /*}}}*/
    13665        int    GeomEdge::Required()       {/*{{{*/
    13766                return type &64;
    13867        }/*}}}*/

  • ../trunk-jpl/src/c/bamg/Mesh.cpp

     
    15931593                        long i0 = GetId(triangles[it][VerticesOfTriangularEdge[j][0]]);
    15941594                        long i1 = GetId(triangles[it][VerticesOfTriangularEdge[j][1]]);
    15951595                        k = edge4->SortAndFind(i0,i1);
    1596                         _assert_(k>=0);
     1596                        if(k<0){
     1597                                delete [] colorV;
     1598                                _error_("This should not happen");
     1599                        }
    15971600                        subdomains[i].direction = (vertices + i0 == edges[k].v[0]) ? 1 : -1;
    15981601                        subdomains[i].edge = edges+k;
    15991602                        Gh.subdomains[i].edge = Gh.edges + k;
    … …  
    20532056                                        if      (xy==0) dd=dxdx;
    20542057                                        else if (xy==1) dd=dxdy;
    20552058                                        else if (xy==2) dd=dydy;
    2056                                         else    _error_("not supported yet");
     2059                                        else{
     2060                                                delete [] detT;
     2061                                                delete [] Mmass;
     2062                                                delete [] workT;
     2063                                                delete [] workV;
     2064                                                delete [] Mmassxx;
     2065                                                delete [] OnBoundary;
     2066                                                _error_("not supported yet");
     2067                                        }
    20572068                                        // do leat 2 iteration for boundary problem
    20582069                                        for (int ijacobi=0;ijacobi<Max(NbJacobi,2);ijacobi++){
    20592070                                                for (i=0;i<nbt;i++)
    … …  
    21022113                delete [] dxdx;
    21032114                delete [] dxdy;
    21042115                delete [] dydy;
    2105                 delete []  workT;
     2116                delete [] workT;
    21062117                delete [] workV;
    21072118                delete [] Mmassxx;
    2108                 delete []  OnBoundary;
     2119                delete [] OnBoundary;
    21092120
    21102121        }
    21112122        /*}}}*/
    … …  
    23932404                        it++;} // end while (it<nbt)
    23942405                        if (nbt == nbtout ||  !NbSubDomTot) {
    23952406                                delete [] HeapArete;
     2407                                delete [] HeapTriangle;
    23962408                                _error_("The boundary is not close: all triangles are outside");
    23972409                        }
    23982410
    … …  
    28822894                                if(kk<10) _printf_("BUG: the geometrical edge " << i << " is on no edge curve\n");
    28832895                        }
    28842896                }
    2885                 if(kk) _error_("See above");
     2897                if(kk){
     2898                        delete [] e;
     2899                        _error_("See above");
     2900                }
    28862901
    28872902                return e;
    28882903        }
    … …  
    43974412                for (i=0;i<Gh.nbv;i++) if (Gh[i].Required()) NbVerticesOnGeomVertex++;
    43984413                printf("\n");
    43994414                if(NbVerticesOnGeomVertex >= maxnbv){
     4415                        delete [] bcurve;
    44004416                        _error_("too many vertices on geometry: " << NbVerticesOnGeomVertex << " >= " << maxnbv);
    44014417                }
    44024418

  • ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

     
    984984        GiaIvinsAnalysisEnum,
    985985        EsaSolutionEnum,
    986986        EsaAnalysisEnum,
    987         MeshdeformationAnalysisEnum,
    988987        LevelsetAnalysisEnum,
    989988        ExtrapolationAnalysisEnum,
    990989        /*}}}*/

  • ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

     
    952952                case GiaIvinsAnalysisEnum : return "GiaIvinsAnalysis";
    953953                case EsaSolutionEnum : return "EsaSolution";
    954954                case EsaAnalysisEnum : return "EsaAnalysis";
    955                 case MeshdeformationAnalysisEnum : return "MeshdeformationAnalysis";
    956955                case LevelsetAnalysisEnum : return "LevelsetAnalysis";
    957956                case ExtrapolationAnalysisEnum : return "ExtrapolationAnalysis";
    958957                case ApproximationEnum : return "Approximation";

  • ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

     
    973973              else if (strcmp(name,"GiaIvinsAnalysis")==0) return GiaIvinsAnalysisEnum;
    974974              else if (strcmp(name,"EsaSolution")==0) return EsaSolutionEnum;
    975975              else if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum;
    976               else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum;
    977976              else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
    978977              else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum;
    979978              else if (strcmp(name,"Approximation")==0) return ApproximationEnum;
    … …  
    997996              else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
    998997              else if (strcmp(name,"MinVel")==0) return MinVelEnum;
    999998              else if (strcmp(name,"MaxVel")==0) return MaxVelEnum;
     999              else if (strcmp(name,"MinVx")==0) return MinVxEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"MinVx")==0) return MinVxEnum;
    1004               else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
     1003              if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
    10051004              else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
    10061005              else if (strcmp(name,"MinVy")==0) return MinVyEnum;
    10071006              else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;

  • ../trunk-jpl/src/c/shared/Numerics/GaussPoints.cpp

     
    16901690        xDelete<IssmPDouble>(work);
    16911691
    16921692}/*}}}*/
    1693 
    1694 /*Element Gauss points TO BE REMOVED*/
    1695 void gaussQuad( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, int nigaus, int njgaus ) {/*{{{*/
    1696         /*Gauss quadrature points for the quadrilaterial.*/
    1697 
    1698         /*get the gauss points using the product of two line rules  */
    1699         GaussLegendreLinear(pxgaus, pxwgt, nigaus);
    1700         GaussLegendreLinear(pegaus, pewgt, njgaus);
    1701 }/*}}}*/
    1702 void gaussHexa( IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, IssmPDouble** pzgaus, IssmPDouble ** pzwgt, int nigaus, int njgaus, int nkgaus ) {/*{{{*/
    1703         /*Gauss quadrature points for the hexahedron.*/
    1704 
    1705         /*  get the gauss points using the product of three line rules  */
    1706         GaussLegendreLinear(pxgaus, pxwgt, nigaus);
    1707         GaussLegendreLinear(pegaus, pewgt, njgaus);
    1708         GaussLegendreLinear(pzgaus, pzwgt, nkgaus);
    1709 }/*}}}*/

  • ../trunk-jpl/src/c/shared/Numerics/GaussPoints.h

     
    1717#define MAX_GAUS_ITER   30
    1818void GaussRecur(IssmPDouble* zero, IssmPDouble* weight, int n, IssmPDouble* alpha, IssmPDouble* beta);
    1919
    20 void gaussQuad(IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, int nigaus, int njgaus);
    21 void gaussHexa(IssmPDouble** pxgaus, IssmPDouble** pxwgt, IssmPDouble** pegaus, IssmPDouble** pewgt, IssmPDouble** pzgaus, IssmPDouble ** pzwgt, int nigaus, int njgaus, int nkgaus);
    22 
    2320#endif

  • ../trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp

     
    6767        }
    6868
    6969}
    70 double ddistance(double x1,double y1,double x2,double y2){
    71         return sqrt(pow(x2-x1,2)+pow(y2-y1,2));
    72 }
    73 double ddot(double x1, double y1, double x2, double y2){
    74         return x1*x2+y1*y2;
    75 }
    76 
    77 bool isPointLeftOfRay(double x, double y, double raySx, double raySy, double rayEx, double rayEy) {
    78           return (y-raySy)*(rayEx-raySx)
    79                     >      (x-raySx)*(rayEy-raySy);
    80 }
    81 
     70double ddistance(double x1,double y1,double x2,double y2){ return sqrt(pow(x2-x1,2)+pow(y2-y1,2)); }
     71double ddot(double x1, double y1, double x2, double y2){ return x1*x2+y1*y2; }
    8272double minimum_distance(double x1, double y1, double x2, double y2, double x0, double y0){
    8373       
    8474        // Return minimum distance between line segment [(x1,y1) (x2,y2)] and point (x0,y0) (v=(x1,y1), w=(x2,y2) and p=(x0,y0)

  • ../trunk-jpl/src/c/Makefile.am

     
    293293                                        ./cores/dummy_core.cpp\
    294294                                        ./cores/surfaceslope_core.cpp\
    295295                                        ./cores/bedslope_core.cpp\
    296                                         ./cores/meshdeformation_core.cpp\
    297296                                        ./cores/damage_core.cpp\
    298297                                        ./cores/levelsetfunctionslope_core.cpp\
    299298                                        ./cores/movingfront_core.cpp\
    … …  
    447446if SMOOTH
    448447issm_sources += ./analyses/SmoothAnalysis.cpp
    449448endif
    450 if MESHDEFORMATION
    451 issm_sources += ./analyses/MeshdeformationAnalysis.cpp
    452 endif
    453449if LEVELSET
    454450issm_sources += ./analyses/LevelsetAnalysis.cpp
    455451issm_sources += ./modules/Calvingx/Calvingx.cpp

  • ../trunk-jpl/src/c/cores/meshdeformation_core.cpp

     
    1 /*!\file: meshdeformation_core.cpp
    2  * \brief: core of the slope solution
    3  */
    4 
    5 #include "./cores.h"
    6 #include "../toolkits/toolkits.h"
    7 #include "../classes/classes.h"
    8 #include "../shared/shared.h"
    9 #include "../solutionsequences/solutionsequences.h"
    10 #include "../modules/modules.h"
    11 
    12 void meshdeformation_core(FemModel* femmodel){
    13 
    14         /*parameters: */
    15         bool save_results;
    16 
    17         /*Recover some parameters: */
    18         femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    19 
    20         if(VerboseSolution()) _printf0_("computing mesh deformation (elasticity)...\n");
    21 
    22         /*Call on core computations: */
    23         femmodel->SetCurrentConfiguration(MeshdeformationAnalysisEnum);
    24         solutionsequence_linear(femmodel);
    25 
    26         if(save_results){
    27                 if(VerboseSolution()) _printf0_("saving results:\n");
    28         }
    29 
    30 }

  • ../trunk-jpl/src/c/classes/Elements/Penta.cpp

     
    12211221        return base*(surface - bed + min( rho_water/rho_ice * bathymetry, 0.) );
    12221222}
    12231223/*}}}*/
    1224 void       Penta::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
    1225 
    1226         /*Intermediary*/
    1227         int    num_controls;
    1228         int*   control_type=NULL;
    1229         Input* input=NULL;
    1230 
    1231         /*retrieve some parameters: */
    1232         this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    1233         this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    1234 
    1235         for(int i=0;i<num_controls;i++){
    1236 
    1237                 if(control_type[i]==MaterialsRheologyBbarEnum){
    1238                         if (!IsOnBase()) goto cleanup_and_return;
    1239                         input=(Input*)this->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input);
    1240                 }
    1241                 else if(control_type[i]==DamageDbarEnum){
    1242                         if (!IsOnBase()) goto cleanup_and_return;
    1243                         input=(Input*)this->inputs->GetInput(DamageDEnum); _assert_(input);
    1244                 }
    1245                 else{
    1246                         input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input);
    1247                 }
    1248                 if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
    1249 
    1250                 ((ControlInput*)input)->UpdateValue(scalar);
    1251                 ((ControlInput*)input)->Constrain();
    1252                 if (save_parameter) ((ControlInput*)input)->SaveValue();
    1253 
    1254                 if(control_type[i]==MaterialsRheologyBbarEnum){
    1255                         this->InputExtrude(MaterialsRheologyBEnum,-1);
    1256                 }
    1257                 else if(control_type[i]==DamageDbarEnum){
    1258                         this->InputExtrude(DamageDEnum,-1);
    1259                 }
    1260         }
    1261 
    1262         /*Clean up and return*/
    1263 cleanup_and_return:
    1264         xDelete<int>(control_type);
    1265 }
    1266 /*}}}*/
    12671224void       Penta::InputDepthAverageAtBase(int original_enum,int average_enum){/*{{{*/
    12681225
    12691226        IssmDouble  Jdet,value;
    … …  
    13671324        }
    13681325}
    13691326/*}}}*/
    1370 void       Penta::InputScale(int enum_type,IssmDouble scale_factor){/*{{{*/
    1371 
    1372         Input* input=NULL;
    1373 
    1374         /*Make a copy of the original input: */
    1375         input=(Input*)this->inputs->GetInput(enum_type);
    1376         if(!input)_error_("could not find old input with enum: " << EnumToStringx(enum_type));
    1377 
    1378         /*Scale: */
    1379         input->Scale(scale_factor);
    1380 }
    1381 /*}}}*/
    13821327void       Penta::InputUpdateFromIoModel(int index,IoModel* iomodel){ /*{{{*/
    13831328
    13841329        /*Intermediaries*/

  • ../trunk-jpl/src/c/classes/Elements/Penta.h

     
    9090                IssmDouble     GroundedArea(void);
    9191                IssmDouble     IceVolume(void);
    9292                IssmDouble     IceVolumeAboveFloatation(void);
    93                 void           InputControlUpdate(IssmDouble scalar,bool save_parameter);
    9493                void           InputDepthAverageAtBase(int enum_type,int average_enum_type);
    9594                void             InputExtrude(int enum_type,int start);
    96                 void           InputScale(int enum_type,IssmDouble scale_factor);
    9795                void           InputUpdateFromIoModel(int index, IoModel* iomodel);
    9896                void           InputUpdateFromSolutionOneDof(IssmDouble* solutiong,int enum_type);
    9997                void           InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solutiong,int enum_type);

  • ../trunk-jpl/src/c/classes/Elements/Seg.h

     
    7878                IssmDouble  GroundedArea(void){_error_("not implemented yet");};
    7979                IssmDouble  IceVolume(void){_error_("not implemented yet");};
    8080                IssmDouble  IceVolumeAboveFloatation(void){_error_("not implemented yet");};
    81                 void        InputControlUpdate(IssmDouble scalar,bool save_parameter){_error_("not implemented yet");};
    8281                void        InputDepthAverageAtBase(int enum_type,int average_enum_type){_error_("not implemented yet");};
    8382                void        InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/};
    84                 void        InputScale(int enum_type,IssmDouble scale_factor){_error_("not implemented yet");};
    8583                void        InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
    8684                void        InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum){_error_("not implemented yet");};
    8785                void        InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum){_error_("not implemented yet");};

  • ../trunk-jpl/src/c/classes/Elements/Tetra.h

     
    9292                bool        IsOnBase();
    9393                bool        IsOnSurface();
    9494                bool        IsZeroLevelset(int levelset_enum){_error_("not implemented");};
    95                 void        InputControlUpdate(IssmDouble scalar,bool save_parameter){_error_("not implemented yet");};
    9695                void        InputDepthAverageAtBase(int enum_type,int average_enum_type){_error_("not implemented yet");};
    9796                void        InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/};
    98                 void        InputScale(int enum_type,IssmDouble scale_factor){_error_("not implemented yet");};
    9997                void        InputUpdateFromIoModel(int index, IoModel* iomodel);
    10098                void        InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum);
    10199                void        InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum){_error_("not implemented yet");};

  • ../trunk-jpl/src/c/classes/Elements/Element.h

     
    221221                virtual IssmDouble GroundedArea(void)=0;
    222222                virtual IssmDouble IceVolume(void)=0;
    223223                virtual IssmDouble IceVolumeAboveFloatation(void)=0;
    224                 virtual void       InputControlUpdate(IssmDouble scalar,bool save_parameter)=0;
    225224                virtual void       InputDepthAverageAtBase(int enum_type,int average_enum_type)=0;
    226225                virtual void       InputExtrude(int input_enum,int start)=0;
    227                 virtual void       InputScale(int enum_type,IssmDouble scale_factor)=0;
    228226                virtual void       InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum)=0;
    229227                virtual void       InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum)=0;
    230228                virtual bool       IsFaceOnBoundary(void)=0;

  • ../trunk-jpl/src/c/classes/Elements/Tria.cpp

     
    13661366        /*Computeportion of the element that has a positive levelset*/
    13671367
    13681368        bool               negative=true;
    1369         int                point;
     1369        int                point=0;
    13701370        const IssmPDouble  epsilon= 1.e-15;
    13711371        IssmDouble         f1,f2;
    13721372
    … …  
    14041404                        f1=gl[1]/(gl[1]-gl[2]);
    14051405                        f2=gl[1]/(gl[1]-gl[0]);
    14061406                }
     1407                else{
     1408                        _error_("This case should NOT be happening");
     1409                }
    14071410        }
    14081411        *point1=point;
    14091412        *fraction1=f1;
    … …  
    16601663        return base*(surface-bed+min(rho_water/rho_ice*bathymetry,0.));
    16611664}
    16621665/*}}}*/
    1663 void       Tria::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
    1664 
    1665         /*Intermediary*/
    1666         int    num_controls;
    1667         int*   control_type=NULL;
    1668         Input* input=NULL;
    1669 
    1670         /*retrieve some parameters: */
    1671         this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    1672         this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    1673 
    1674         for(int i=0;i<num_controls;i++){
    1675                 input=(Input*)this->inputs->GetInput(control_type[i]);   _assert_(input);
    1676                 if (input->ObjectEnum()!=ControlInputEnum){
    1677                         _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
    1678                 }
    1679 
    1680                 ((ControlInput*)input)->UpdateValue(scalar);
    1681                 ((ControlInput*)input)->Constrain();
    1682                 if (save_parameter) ((ControlInput*)input)->SaveValue();
    1683 
    1684         }
    1685 
    1686         /*Clean up and return*/
    1687         xDelete<int>(control_type);
    1688 }
    1689 /*}}}*/
    16901666void       Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){/*{{{*/
    16911667
    16921668        /*New input*/
    … …  
    17051681        this->inputs->AddInput((Input*)newinput);
    17061682}
    17071683/*}}}*/
    1708 void       Tria::InputScale(int enum_type,IssmDouble scale_factor){/*{{{*/
    1709 
    1710         Input* input=NULL;
    1711 
    1712         /*Make a copy of the original input: */
    1713         input=(Input*)this->inputs->GetInput(enum_type);
    1714         if(!input)_error_("could not find old input with enum: " << EnumToStringx(enum_type));
    1715 
    1716         /*Scale: */
    1717         input->Scale(scale_factor);
    1718 }
    1719 /*}}}*/
    17201684void       Tria::InputUpdateFromIoModel(int index, IoModel* iomodel){ //i is the element index/*{{{*/
    17211685
    17221686        /*Intermediaries*/

  • ../trunk-jpl/src/c/classes/Elements/Tria.h

     
    9090                bool        HasEdgeOnSurface();
    9191                IssmDouble  IceVolume(void);
    9292                IssmDouble  IceVolumeAboveFloatation(void);
    93                 void        InputControlUpdate(IssmDouble scalar,bool save_parameter);
    9493                void        InputDepthAverageAtBase(int enum_type,int average_enum_type);
    9594                void        InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/};
    96                 void        InputScale(int enum_type,IssmDouble scale_factor);
    9795                bool            IsFaceOnBoundary(void);
    9896                bool            IsIcefront(void);
    9997                bool        IsNodeOnShelfFromFlags(IssmDouble* flags);

  • ../trunk-jpl/src/c/classes/Inputs/ControlInput.h

     
    8080                void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");};
    8181                void SaveValue(void);
    8282                void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
    83                 void ScaleGradient(IssmDouble scale);
    8483                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    8584                void SetGradient(Input* gradient_in);
    8685                void SetInput(Input* in_input);

  • ../trunk-jpl/src/c/classes/Inputs/DatasetInput.h

     
    7676                void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");};
    7777                void SaveValue(void){_error_("not implemented yet");};
    7878                void Scale(IssmDouble scale_factor){_error_("not implemented yet");};
    79                 void ScaleGradient(IssmDouble scale){_error_("not implemented yet");};
    8079                void Set(IssmDouble setvalue){_error_("Set not implemented yet");};
    8180                void SetGradient(Input* gradient_in){_error_("not implemented yet");};
    8281                void UpdateValue(IssmDouble scalar){_error_("not implemented yet");};

  • ../trunk-jpl/src/c/classes/Inputs/ControlInput.cpp

     
    235235        if(savedvalues) delete this->savedvalues;
    236236        this->savedvalues=xDynamicCast<Input*>(this->values->copy());
    237237}/*}}}*/
    238 void ControlInput::ScaleGradient(IssmDouble scaling_factor){/*{{{*/
    239         if(!gradient) _error_("Gradient of ControlInput " << EnumToStringx(enum_type) << " not found");
    240         gradient->Scale(scaling_factor);
    241 }/*}}}*/
    242238void ControlInput::SetGradient(Input* gradient_in){/*{{{*/
    243239
    244240        /*Get enum for current gradient*/

  • ../trunk-jpl/src/c/analyses/MeshdeformationAnalysis.cpp

     
    1 #include "./MeshdeformationAnalysis.h"
    2 #include "../toolkits/toolkits.h"
    3 #include "../classes/classes.h"
    4 #include "../shared/shared.h"
    5 #include "../modules/modules.h"
    6 
    7 /*Model processing*/
    8 void MeshdeformationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
    9            _error_("not implemented yet");
    10 }/*}}}*/
    11 void MeshdeformationAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
    12            _error_("not implemented yet");
    13 }/*}}}*/
    14 void MeshdeformationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
    15            _error_("not implemented yet");
    16 }/*}}}*/
    17 int  MeshdeformationAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
    18         _error_("not implemented");
    19 }/*}}}*/
    20 void MeshdeformationAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    21            _error_("not implemented yet");
    22 }/*}}}*/
    23 void MeshdeformationAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    24            _error_("not implemented yet");
    25 }/*}}}*/
    26 
    27 /*Finite Element Analysis*/
    28 void           MeshdeformationAnalysis::Core(FemModel* femmodel){/*{{{*/
    29         _error_("not implemented");
    30 }/*}}}*/
    31 ElementVector* MeshdeformationAnalysis::CreateDVector(Element* element){/*{{{*/
    32         /*Default, return NULL*/
    33         return NULL;
    34 }/*}}}*/
    35 ElementMatrix* MeshdeformationAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
    36 _error_("Not implemented");
    37 }/*}}}*/
    38 ElementMatrix* MeshdeformationAnalysis::CreateKMatrix(Element* element){/*{{{*/
    39         _error_("not implemented yet");
    40 }/*}}}*/
    41 ElementVector* MeshdeformationAnalysis::CreatePVector(Element* element){/*{{{*/
    42 _error_("not implemented yet");
    43 }/*}}}*/
    44 void           MeshdeformationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    45            _error_("not implemented yet");
    46 }/*}}}*/
    47 void           MeshdeformationAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
    48         _error_("Not implemented yet");
    49 }/*}}}*/
    50 void           MeshdeformationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    51         _error_("not implemented yet");
    52 }/*}}}*/
    53 void           MeshdeformationAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
    54         /*Default, do nothing*/
    55         return;
    56 }/*}}}*/

  • ../trunk-jpl/src/c/analyses/MeshdeformationAnalysis.h

     
    1 /*! \file MeshdeformationAnalysis.h
    2  *  \brief: header file for generic external result object
    3  */
    4 
    5 #ifndef _MeshdeformationAnalysis_
    6 #define _MeshdeformationAnalysis_
    7 
    8 /*Headers*/
    9 #include "./Analysis.h"
    10 
    11 class MeshdeformationAnalysis: public Analysis{
    12 
    13         public:
    14                 /*Model processing*/
    15                 void CreateConstraints(Constraints* constraints,IoModel* iomodel);
    16                 void CreateLoads(Loads* loads, IoModel* iomodel);
    17                 void CreateNodes(Nodes* nodes,IoModel* iomodel);
    18                 int  DofsPerNode(int** doflist,int domaintype,int approximation);
    19                 void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
    20                 void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
    21 
    22                 /*Finite element Analysis*/
    23                 void           Core(FemModel* femmodel);
    24                 ElementVector* CreateDVector(Element* element);
    25                 ElementMatrix* CreateJacobianMatrix(Element* element);
    26                 ElementMatrix* CreateKMatrix(Element* element);
    27                 ElementVector* CreatePVector(Element* element);
    28                 void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
    30                 void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    31                 void           UpdateConstraints(FemModel* femmodel);
    32 };
    33 #endif

  • ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

     
    551551                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    552552        }
    553553}/*}}}*/
    554 void           LevelsetAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
    555 
    556         /* Intermediaries */
    557         int i,k;
    558         IssmDouble dmaxp=0.,dmaxm=0,val=0.;
    559 
    560         /*Initialize vector with number of vertices*/
    561         int numvertices=femmodel->vertices->NumberOfVertices();
    562         Element* element;
    563 
    564         Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
    565         GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexPIdEnum);
    566        
    567         /* set NaN on elements intersected by zero levelset */
    568         for(i=0;i<femmodel->elements->Size();i++){
    569                 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    570                 if(element->IsZeroLevelset(MaskIceLevelsetEnum))
    571                         for(k=0;k<element->GetNumberOfVertices();k++)
    572                                 vec_dist_zerolevelset->SetValue(element->vertices[k]->Sid(),NAN,INS_VAL);
    573         }
    574 
    575         /* set distance on elements intersected by zero levelset */
    576         for(i=0;i<femmodel->elements->Size();i++){
    577                 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    578                 if(element->IsZeroLevelset(MaskIceLevelsetEnum)){
    579                         SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
    580 
    581                         /* Get maximum distance to interface along vertices */
    582                         for(k=0;k<element->GetNumberOfVertices();k++){
    583                                         vec_dist_zerolevelset->GetValue(&val,element->vertices[k]->Sid());
    584                                         if((val>0.) && (val>dmaxp))
    585                                                  dmaxp=val;
    586                                         else if((val<0.) && (val<dmaxm))
    587                                                  dmaxm=val;
    588                         }
    589                 }
    590         }
    591 
    592         /* set all none intersected vertices to max/min distance */
    593         for(i=0;i<numvertices;i++){
    594                 vec_dist_zerolevelset->GetValue(&val,i);
    595                 if(val==1.) //FIXME: improve check
    596                         vec_dist_zerolevelset->SetValue(i,3.*dmaxp,INS_VAL);
    597                 else if(val==-1.)
    598                         vec_dist_zerolevelset->SetValue(i,3.*dmaxm,INS_VAL);
    599         }
    600 
    601         /*Assemble vector and serialize */
    602         vec_dist_zerolevelset->Assemble();
    603         IssmDouble* dist_zerolevelset=vec_dist_zerolevelset->ToMPISerial();
    604         InputUpdateFromVectorx(femmodel,dist_zerolevelset,MaskIceLevelsetEnum,VertexSIdEnum);
    605 
    606         /*Clean up and return*/
    607         delete vec_dist_zerolevelset;
    608         delete dist_zerolevelset;
    609 }/*}}}*/
    610554void           LevelsetAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
    611555
    612556        if(!element->IsZeroLevelset(MaskIceLevelsetEnum))

  • ../trunk-jpl/src/c/analyses/analyses.h

     
    3333#include "./SmbAnalysis.h"
    3434#include "./SealevelriseAnalysis.h"
    3535#include "./MeltingAnalysis.h"
    36 #include "./MeshdeformationAnalysis.h"
    3736#include "./SmoothAnalysis.h"
    3837#include "./StressbalanceAnalysis.h"
    3938#include "./StressbalanceSIAAnalysis.h"

  • ../trunk-jpl/src/c/analyses/LevelsetAnalysis.h

     
    3131                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    3232                void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
    3333                void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
    34                 void           SetDistanceOnIntersectedElements(FemModel* femmodel);
    3534                void           SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element);
    3635                void           UpdateConstraints(FemModel* femmodel);
    3736};

  • ../trunk-jpl/src/c/analyses/EnumToAnalysis.cpp

     
    109109                #ifdef _HAVE_ESA_
    110110                case EsaAnalysisEnum : return new EsaAnalysis();
    111111                #endif
    112                 #ifdef _HAVE_MESHDEFORMATION_
    113                 case MeshdeformationAnalysisEnum : return new MeshdeformationAnalysis();
    114                 #endif
    115112                #ifdef _HAVE_LEVELSET_
    116113                case LevelsetAnalysisEnum : return new LevelsetAnalysis();
    117114                #endif
Note: See TracBrowser for help on using the repository browser.