Changeset 5635


Ignore:
Timestamp:
08/31/10 12:15:10 (15 years ago)
Author:
Mathieu Morlighem
Message:

more GaussTria

Location:
issm/trunk/src/c/objects/Elements
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r5633 r5635  
    3232                virtual void   GetSolutionFromInputs(Vec solution)=0;
    3333                virtual void   GetNodes(void** nodes)=0;
    34                 virtual void*  GetMatPar()=0;
    3534                virtual bool   GetShelf()=0;
    3635                virtual bool   GetOnBed()=0;

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

    r5633 r5635  
    941941        this->results=new Results();
    942942
    943 }
    944 /*}}}*/
    945 /*FUNCTION Penta::GetMatPar {{{1*/
    946 void* Penta::GetMatPar(){
    947 
    948         return matpar;
    949943}
    950944/*}}}*/
    … …  
    48604854
    48614855        /*clean-up*/
    4862         delete gauss;
     4856        //delete gauss;
    48634857}
    48644858/*}}}*/

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

    r5633 r5635  
    7979                void   CreatePVector(Vec pg);
    8080                void   DeleteResults(void);
    81                 void*  GetMatPar();
    8281                void   GetNodes(void** nodes);
    8382                bool   GetOnBed();

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

    r5634 r5635  
    572572/*FUNCTION Tria::ComputeBasalStress {{{1*/
    573573void  Tria::ComputeBasalStress(Vec eps){
    574 
    575         int i;
    576         const int numvertices=3;
    577         int doflist[numvertices];
    578         double value;
    579 
    580         /*plug local pressure values into global pressure vector: */
    581574        ISSMERROR("Not Implemented yet");
    582         //VecSetValues(eps,1,2,(const double*)&value,INSERT_VALUES);
    583 
    584575}
    585576/*}}}*/
    … …  
    587578void  Tria::ComputePressure(Vec pg){
    588579
    589         int i;
    590         const int numvertices=3;
    591         int doflist[numvertices];
    592         double pressure[numvertices];
    593         double thickness[numvertices];
    594         double rho_ice,g;
     580        const int numvertices= 3;
     581        int       doflist[numvertices];
     582        double    pressure[numvertices];
     583        double    thickness[numvertices];
     584        double    rho_ice,g;
    595585       
    596586        /*Get dof list on which we will plug the pressure values: */
    … …  
    602592        GetParameterListOnVertices(&thickness[0],ThicknessEnum);
    603593
    604         for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
     594        for(int i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
    605595
    606596        /*plug local pressure values into global pressure vector: */
    … …  
    611601/*FUNCTION Tria::ComputeStrainRate {{{1*/
    612602void  Tria::ComputeStrainRate(Vec eps){
    613 
    614         int       i;
    615         const int numvertices = 3;
    616         int       doflist[numvertices];
    617         double    value;
    618 
    619         /*plug local pressure values into global pressure vector: */
    620603        ISSMERROR("Not Implemented yet");
    621         //VecSetValues(eps,1,2,(const double*)&value,INSERT_VALUES);
    622 
    623604}
    624605/*}}}*/
    … …  
    843824}
    844825/*}}}*/
    845 /*FUNCTION Tria::GetMatPar {{{1*/
    846 void* Tria::GetMatPar(){
    847 
    848         return matpar;
    849 }
    850 /*}}}*/
    851826/*FUNCTION Tria::GetNodes {{{1*/
    852827void  Tria::GetNodes(void** vpnodes){
    … …  
    13071282        }
    13081283
    1309         /*clean up*/
     1284        /*clean up and return*/
    13101285        xfree((void**)&new_inputs);
    13111286        xfree((void**)&old_inputs);
    1312 
    1313         /*Return output*/
    13141287        return converged;
    13151288
    … …  
    17831756double Tria::ThicknessAbsMisfit(bool process_units){
    17841757
    1785         int i;
    1786 
    1787         /* output: */
    1788         double Jelem=0;
    1789 
    1790         /* node data: */
     1758        /* Constants */
    17911759        const int    numvertices=3;
    17921760        const int    numdof=1*numvertices;
    1793         const int    NDOF=1;
    1794         double       xyz_list[numvertices][3];
    1795 
    1796         /* gaussian points: */
    1797         int     num_gauss,ig;
    1798         double* first_gauss_area_coord  =  NULL;
    1799         double* second_gauss_area_coord =  NULL;
    1800         double* third_gauss_area_coord  =  NULL;
    1801         double* gauss_weights           =  NULL;
    1802         double  gauss_weight;
    1803         double  gauss_l1l2l3[3];
    1804 
    1805         /* parameters: */
    1806         double  thickness,thicknessobs,weight;
    1807 
    1808         /* Jacobian: */
    1809         double Jdet;
    1810 
    1811         /*inputs: */
    1812         bool onwater;
    1813         Input* thickness_input=NULL;
    1814         Input* thicknessobs_input=NULL;
    1815         Input* weights_input=NULL;
    1816 
    1817         /*retrieve inputs :*/
    1818         inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1761
     1762        /*Intermediaries*/
     1763        int        i,ig;
     1764        bool       onwater;
     1765        double     thickness,thicknessobs,weight;
     1766        double     Jdet;
     1767        double     Jelem = 0;
     1768        double     xyz_list[numvertices][3];
     1769        GaussTria *gauss = NULL;
    18191770
    18201771        /*If on water, return 0: */
    … …  
    18251776
    18261777        /*Retrieve all inputs we will be needing: */
    1827         thickness_input   =inputs->GetInput(ThicknessEnum);   ISSMASSERT(thickness_input);
    1828         thicknessobs_input=inputs->GetInput(ThicknessObsEnum);ISSMASSERT(thicknessobs_input);
    1829         weights_input     =inputs->GetInput(WeightsEnum);ISSMASSERT(weights_input);
    1830 
    1831         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    1832         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     1778        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     1779        Input* thickness_input   =inputs->GetInput(ThicknessEnum);   ISSMASSERT(thickness_input);
     1780        Input* thicknessobs_input=inputs->GetInput(ThicknessObsEnum);ISSMASSERT(thicknessobs_input);
     1781        Input* weights_input     =inputs->GetInput(WeightsEnum);ISSMASSERT(weights_input);
    18331782
    18341783        /* Start  looping on the number of gaussian points: */
    1835         for (ig=0; ig<num_gauss; ig++){
    1836                 /*Pick up the gaussian point: */
    1837                 gauss_weight=*(gauss_weights+ig);
    1838                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    1839                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    1840                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     1784        gauss=new GaussTria(2);
     1785        for (ig=gauss->begin();ig<gauss->end();ig++){
     1786
     1787                gauss->GaussPoint(ig);
    18411788
    18421789                /* Get Jacobian determinant: */
    1843                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     1790                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    18441791
    18451792                /*Get parameters at gauss point*/
    1846                 thickness_input->GetParameterValue(&thickness, &gauss_l1l2l3[0]);
    1847                 thicknessobs_input->GetParameterValue(&thicknessobs, &gauss_l1l2l3[0]);
    1848                 weights_input->GetParameterValue(&weight, &gauss_l1l2l3[0]);
     1793                thickness_input->GetParameterValue(&thickness,gauss);
     1794                thicknessobs_input->GetParameterValue(&thicknessobs,gauss);
     1795                weights_input->GetParameterValue(&weight,gauss);
    18491796
    18501797                /*compute ThicknessAbsMisfit*/
    1851                 Jelem+=0.5*pow(thickness-thicknessobs,2.0)*weight*Jdet*gauss_weight;
    1852         }
    1853 
    1854         xfree((void**)&first_gauss_area_coord);
    1855         xfree((void**)&second_gauss_area_coord);
    1856         xfree((void**)&third_gauss_area_coord);
    1857         xfree((void**)&gauss_weights);
    1858 
    1859         /*Return: */
     1798                Jelem+=0.5*pow(thickness-thicknessobs,2.0)*weight*Jdet*gauss->weight;
     1799        }
     1800
     1801        /* clean up and Return: */
     1802        delete gauss;
    18601803        return Jelem;
    18611804}
    … …  
    18851828
    18861829        /* gaussian points: */
    1887         int     num_gauss,ig;
    1888         double* first_gauss_area_coord  =  NULL;
    1889         double* second_gauss_area_coord =  NULL;
    1890         double* third_gauss_area_coord  =  NULL;
    1891         double* gauss_weights           =  NULL;
    1892         double  gauss_weight;
    1893         double  gauss_l1l2l3[3];
     1830        int     ig;
     1831        GaussTria *gauss=NULL;
    18941832
    18951833        /* parameters: */
    … …  
    19531891        }
    19541892
    1955         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    1956         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    1957 
    19581893        /* Start  looping on the number of gaussian points: */
    1959         for (ig=0; ig<num_gauss; ig++){
    1960                 /*Pick up the gaussian point: */
    1961                 gauss_weight=*(gauss_weights+ig);
    1962                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    1963                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    1964                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    1965 
    1966                 /* Get Jacobian determinant: */
    1967                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    1968 
    1969                 /*Compute misfit at gaussian point: */
    1970                 TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss_l1l2l3);
    1971 
    1972                 /*compute SurfaceAbsVelMisfit*/
    1973                 Jelem+=misfit*Jdet*gauss_weight;
    1974         }
    1975 
    1976         xfree((void**)&first_gauss_area_coord);
    1977         xfree((void**)&second_gauss_area_coord);
    1978         xfree((void**)&third_gauss_area_coord);
    1979         xfree((void**)&gauss_weights);
    1980                
    1981         /*Return: */
     1894        gauss=new GaussTria(2);
     1895        for (ig=gauss->begin();ig<gauss->end(); ig++){
     1896
     1897                gauss->GaussPoint(ig);
     1898
     1899                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     1900                TriaRef::GetParameterValue(&misfit, &misfit_list[0],gauss);
     1901                Jelem+=misfit*Jdet*gauss->weight;
     1902        }
     1903
     1904        /*clean up and Return: */
     1905        delete gauss;
    19821906        return Jelem;
    19831907}

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

    r5633 r5635  
    7575                void   CreateKMatrix(Mat Kgg);
    7676                void   CreatePVector(Vec pg);
    77                 void*  GetMatPar();
    7877                void   GetNodes(void** nodes);
    7978                bool   GetOnBed();
Note: See TracChangeset for help on using the changeset viewer.