Changeset 5286


Ignore:
Timestamp:
08/16/10 13:28:18 (15 years ago)
Author:
Mathieu Morlighem
Message:

Fixed enormous bug in Misfitx and added c/modules/ThicknessAbsMisfitx needed for Balancedthicknesses"

Location:
issm/trunk/src/c
Files:
3 added
11 edited

Legend:

Unmodified
Added
Removed

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

    r5281 r5286  
    198198        KflagEnum,
    199199        MassFluxEnum,
     200        ThicknessAbsMisfitEnum,
    200201        SurfaceAbsVelMisfitEnum,
    201202        SurfaceRelVelMisfitEnum,

  • issm/trunk/src/c/EnumDefinitions/EnumToString.cpp

    r5281 r5286  
    172172                case KflagEnum : return "Kflag";
    173173                case MassFluxEnum : return "MassFlux";
     174                case ThicknessAbsMisfitEnum : return "ThicknessAbsMisfit";
    174175                case SurfaceAbsVelMisfitEnum : return "SurfaceAbsVelMisfit";
    175176                case SurfaceRelVelMisfitEnum : return "SurfaceRelVelMisfit";

  • issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp

    r5283 r5286  
    170170        else if (strcmp(name,"Kflag")==0) return KflagEnum;
    171171        else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
     172        else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
    172173        else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
    173174        else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;

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

    r5281 r5286  
    443443                                        ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h\
    444444                                        ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp\
     445                                        ./modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.h\
     446                                        ./modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.cpp\
    445447                                        ./modules/CostFunctionx/CostFunctionx.h\
    446448                                        ./modules/CostFunctionx/CostFunctionx.cpp\
    … …  
    986988                                        ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.h\
    987989                                        ./modules/SurfaceAverageVelMisfitx/SurfaceAverageVelMisfitx.cpp\
     990                                        ./modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.h\
     991                                        ./modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.cpp\
    988992                                        ./modules/CostFunctionx/CostFunctionx.h\
    989993                                        ./modules/CostFunctionx/CostFunctionx.cpp\

  • issm/trunk/src/c/modules/Responsex/Responsex.cpp

    r5281 r5286  
    5151                MaxAbsVzx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units);
    5252        }
     53        else if(strncmp(response_descriptor,"MassFlux",8)==0){
     54                MassFluxx(presponse,elements,nodes,vertices,loads,materials,parameters,response_descriptor,process_units);
     55        }
    5356        else if(strcmp(response_descriptor,"SurfaceAbsVelMisfit")==0){
    5457                SurfaceAbsVelMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units);
    … …  
    6669                SurfaceAverageVelMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units);
    6770        }
    68         else if(strncmp(response_descriptor,"MassFlux",8)==0){
    69                 MassFluxx(presponse,elements,nodes,vertices,loads,materials,parameters,response_descriptor,process_units);
     71        else if(strcmp(response_descriptor,"ThicknessAbsMisfit")==0){
     72                ThicknessAbsMisfitx( presponse, elements,nodes, vertices, loads, materials, parameters,process_units);
    7073        }
    7174        else{

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

    r5281 r5286  
    5555#include "./MinVyx/MinVyx.h"
    5656#include "./MinVzx/MinVzx.h"
     57#include "./ThicknessAbsMisfitx/ThicknessAbsMisfitx.h"
    5758#include "./SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h"
    5859#include "./SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h"

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

    r5281 r5286  
    4242                virtual void   GradjDrag(Vec gradient)=0;
    4343                virtual void   GradjB(Vec gradient)=0;
     44                virtual double ThicknessAbsMisfit(bool process_units)=0;
    4445                virtual double SurfaceAbsVelMisfit(bool process_units)=0;
    4546                virtual double SurfaceRelVelMisfit(bool process_units)=0;

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

    r5281 r5286  
    15971597}
    15981598/*}}}*/
    1599 /*FUNCTION Penta::SurfaceAbsVelMisfit {{{1*/
    1600 double Penta::SurfaceAbsVelMisfit(bool process_units){
     1599/*FUNCTION Penta::ThicknessAbsMisfit {{{1*/
     1600double Penta::ThicknessAbsMisfit(bool process_units){
    16011601
    16021602        double J;
    … …  
    16171617        /*If on water, return 0: */
    16181618        if(onwater)return 0;
     1619        ISSMERROR("Not implemented yet");
     1620
     1621        tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
     1622        J=tria->ThicknessAbsMisfit(process_units);
     1623        delete tria->matice; delete tria;
     1624        return J;
     1625}
     1626/*}}}*/
     1627/*FUNCTION Penta::SurfaceAbsVelMisfit {{{1*/
     1628double Penta::SurfaceAbsVelMisfit(bool process_units){
     1629
     1630        double J;
     1631        Tria* tria=NULL;
     1632
     1633        /*inputs: */
     1634        bool onwater;
     1635        bool onsurface;
     1636        bool onbed;
     1637        int  approximation;
     1638
     1639        /*retrieve inputs :*/
     1640        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1641        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     1642        inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
     1643        inputs->GetParameterValue(&approximation,ApproximationEnum);
     1644
     1645        /*If on water, return 0: */
     1646        if(onwater)return 0;
    16191647
    16201648        /*Bail out if this element if:
    16211649         * -> Non MacAyeal and not on the surface
    16221650         * -> MacAyeal (2d model) and not on bed) */
    1623         if ((approximation!=MacAyealApproximationEnum && !onsurface) || (MacAyealApproximationEnum && !onbed)){
     1651        if ((approximation!=MacAyealApproximationEnum && !onsurface) || (approximation==MacAyealApproximationEnum && !onbed)){
    16241652                return 0;
    16251653        }
    1626         else if (MacAyealApproximationEnum){
     1654        else if (approximation==MacAyealApproximationEnum){
    16271655
    16281656                /*This element should be collapsed into a tria element at its base. Create this tria element,

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

    r5281 r5286  
    104104                void   MinVy(double* pminvy, bool process_units);
    105105                void   MinVz(double* pminvz, bool process_units);
     106                double ThicknessAbsMisfit(bool process_units);
    106107                double SurfaceAbsVelMisfit(bool process_units);
    107108                double SurfaceRelVelMisfit(bool process_units);

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

    r5281 r5286  
    573573        double  gauss_weight;
    574574        double  gauss_l1l2l3[3];
    575         double  k_gauss;
    576575        double  B_gauss;
     576        double  dhdt_gauss;
    577577
    578578        /* parameters: */
    … …  
    627627                else if (control_type==RheologyB2dEnum){
    628628                        matice->inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyB2dEnum);
     629                        Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
     630                }
     631                else if (control_type==DhDtEnum){
     632                        matice->inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],DhDtEnum);
    629633                        Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
    630634                }
    … …  
    18551859        *pminvz=minvz;
    18561860
     1861}
     1862/*}}}*/
     1863/*FUNCTION Tria::ThicknessAbsMisfit {{{1*/
     1864double Tria::ThicknessAbsMisfit(bool process_units){
     1865
     1866        int i;
     1867
     1868        /* output: */
     1869        double Jelem=0;
     1870
     1871        /* node data: */
     1872        const int    numgrids=3;
     1873        const int    numdof=1*numgrids;
     1874        const int    NDOF=1;
     1875        double       xyz_list[numgrids][3];
     1876
     1877        /* gaussian points: */
     1878        int     num_gauss,ig;
     1879        double* first_gauss_area_coord  =  NULL;
     1880        double* second_gauss_area_coord =  NULL;
     1881        double* third_gauss_area_coord  =  NULL;
     1882        double* gauss_weights           =  NULL;
     1883        double  gauss_weight;
     1884        double  gauss_l1l2l3[3];
     1885        double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     1886
     1887        /* parameters: */
     1888        double  thickness,thicknessobs,weight;
     1889
     1890        /* Jacobian: */
     1891        double Jdet;
     1892
     1893        /*inputs: */
     1894        bool onwater;
     1895        Input* thickness_input=NULL;
     1896        Input* thicknessobs_input=NULL;
     1897        Input* weights_input=NULL;
     1898
     1899        /*retrieve inputs :*/
     1900        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1901
     1902        /*If on water, return 0: */
     1903        if(onwater)return 0;
     1904
     1905        /* Get node coordinates and dof list: */
     1906        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     1907
     1908        /*Retrieve all inputs we will be needing: */
     1909        thickness_input   =inputs->GetInput(ThicknessEnum);   ISSMASSERT(thickness_input);
     1910        thicknessobs_input=inputs->GetInput(ThicknessObsEnum);ISSMASSERT(thicknessobs_input);
     1911        weights_input     =inputs->GetInput(WeightsEnum);ISSMASSERT(weights_input);
     1912
     1913        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     1914        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     1915
     1916        /* Start  looping on the number of gaussian points: */
     1917        for (ig=0; ig<num_gauss; ig++){
     1918                /*Pick up the gaussian point: */
     1919                gauss_weight=*(gauss_weights+ig);
     1920                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     1921                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     1922                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     1923
     1924                /* Get Jacobian determinant: */
     1925                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     1926
     1927                /*Get parameters at gauss point*/
     1928                thickness_input->GetParameterValue(&thickness, &gauss_l1l2l3[0]);
     1929                thicknessobs_input->GetParameterValue(&thicknessobs, &gauss_l1l2l3[0]);
     1930                weights_input->GetParameterValue(&weight, &gauss_l1l2l3[0]);
     1931
     1932                /*compute ThicknessAbsMisfit*/
     1933                Jelem+=0.5*pow(thickness_input-thicknessobs_input,2.0)*weight*Jdet*gauss_weight;
     1934        }
     1935
     1936        xfree((void**)&first_gauss_area_coord);
     1937        xfree((void**)&second_gauss_area_coord);
     1938        xfree((void**)&third_gauss_area_coord);
     1939        xfree((void**)&gauss_weights);
     1940
     1941        /*Return: */
     1942        return Jelem;
    18571943}
    18581944/*}}}*/
    … …  
    48054891        thickness_input   =inputs->GetInput(ThicknessEnum);   ISSMASSERT(thickness_input);
    48064892        thicknessobs_input=inputs->GetInput(ThicknessObsEnum);ISSMASSERT(thicknessobs_input);
    4807         weights_input     =inputs->GetInput(ThicknessObsEnum);ISSMASSERT(weights_input);
     4893        weights_input     =inputs->GetInput(WeightsEnum);ISSMASSERT(weights_input);
    48084894
    48094895        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */

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

    r5281 r5286  
    103103                void   MinVy(double* pminvy, bool process_units);
    104104                void   MinVz(double* pminvz, bool process_units);
     105                double ThicknessAbsMisfit(bool process_units);
    105106                double SurfaceAbsVelMisfit(bool process_units);
    106107                double SurfaceRelVelMisfit(bool process_units);
Note: See TracChangeset for help on using the changeset viewer.