Ignore:
Timestamp:
09/16/11 08:23:21 (14 years ago)
Author:
Eric.Larour
Message:

Conditional response functionality from elements

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r9777 r9817  
    19261926}
    19271927/*}}}*/
    1928 /*FUNCTION Tria::MassFlux {{{1*/
    1929 double Tria::MassFlux( double* segment,bool process_units){
    1930 
    1931         const int    numdofs=2;
    1932 
    1933         int        i;
    1934         double     mass_flux=0;
    1935         double     xyz_list[NUMVERTICES][3];
    1936         double     normal[2];
    1937         double     length,rho_ice;
    1938         double     x1,y1,x2,y2,h1,h2;
    1939         double     vx1,vx2,vy1,vy2;
    1940         GaussTria* gauss_1=NULL;
    1941         GaussTria* gauss_2=NULL;
    1942 
    1943         /*Get material parameters :*/
    1944         rho_ice=matpar->GetRhoIce();
    1945 
    1946         /*First off, check that this segment belongs to this element: */
    1947         if ((int)*(segment+4)!=this->id)_error_("%s%i%s%i","error message: segment with id ",(int)*(segment+4)," does not belong to element with id:",this->id);
    1948 
    1949         /*Recover segment node locations: */
    1950         x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
    1951 
    1952         /*Get xyz list: */
    1953         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    1954 
    1955         /*get area coordinates of 0 and 1 locations: */
    1956         gauss_1=new GaussTria();
    1957         gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
    1958         gauss_2=new GaussTria();
    1959         gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
    1960 
    1961         normal[0]=cos(atan2(x1-x2,y2-y1));
    1962         normal[1]=sin(atan2(x1-x2,y2-y1));
    1963 
    1964         length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2));
    1965 
    1966         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    1967         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    1968         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    1969 
    1970         thickness_input->GetParameterValue(&h1, gauss_1);
    1971         thickness_input->GetParameterValue(&h2, gauss_2);
    1972         vx_input->GetParameterValue(&vx1,gauss_1);
    1973         vx_input->GetParameterValue(&vx2,gauss_2);
    1974         vy_input->GetParameterValue(&vy1,gauss_1);
    1975         vy_input->GetParameterValue(&vy2,gauss_2);
    1976 
    1977         mass_flux= rho_ice*length*( 
    1978                                 (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
    1979                                 (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
    1980                                 );
    1981 
    1982         /*Process units: */
    1983         mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum);
    1984 
    1985         /*clean up and return:*/
    1986         delete gauss_1;
    1987         delete gauss_2;
    1988         return mass_flux;
    1989 }
    1990 /*}}}*/
    1991 /*FUNCTION Tria::MaxAbsVx{{{1*/
    1992 void  Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){
    1993 
    1994         /*Get maximum:*/
    1995         double maxabsvx=this->inputs->MaxAbs(VxEnum);
    1996 
    1997         /*process units if requested: */
    1998         if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum);
    1999 
    2000         /*Assign output pointers:*/
    2001         *pmaxabsvx=maxabsvx;
    2002 }
    2003 /*}}}*/
    2004 /*FUNCTION Tria::MaxAbsVy{{{1*/
    2005 void  Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){
    2006 
    2007         /*Get maximum:*/
    2008         double maxabsvy=this->inputs->MaxAbs(VyEnum);
    2009 
    2010         /*process units if requested: */
    2011         if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum);
    2012 
    2013         /*Assign output pointers:*/
    2014         *pmaxabsvy=maxabsvy;
    2015 }
    2016 /*}}}*/
    2017 /*FUNCTION Tria::MaxAbsVz{{{1*/
    2018 void  Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){
    2019 
    2020         /*Get maximum:*/
    2021         double maxabsvz=this->inputs->MaxAbs(VzEnum);
    2022 
    2023         /*process units if requested: */
    2024         if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum);
    2025 
    2026         /*Assign output pointers:*/
    2027         *pmaxabsvz=maxabsvz;
    2028 }
    2029 /*}}}*/
    2030 /*FUNCTION Tria::MaxVel{{{1*/
    2031 void  Tria::MaxVel(double* pmaxvel, bool process_units){
    2032 
    2033         /*Get maximum:*/
    2034         double maxvel=this->inputs->Max(VelEnum);
    2035 
    2036         /*process units if requested: */
    2037         if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum);
    2038 
    2039         /*Assign output pointers:*/
    2040         *pmaxvel=maxvel;
    2041 }
    2042 /*}}}*/
    2043 /*FUNCTION Tria::MaxVx{{{1*/
    2044 void  Tria::MaxVx(double* pmaxvx, bool process_units){
    2045 
    2046         /*Get maximum:*/
    2047         double maxvx=this->inputs->Max(VxEnum);
    2048 
    2049         /*process units if requested: */
    2050         if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum);
    2051 
    2052         /*Assign output pointers:*/
    2053         *pmaxvx=maxvx;
    2054 }
    2055 /*}}}*/
    2056 /*FUNCTION Tria::MaxVy{{{1*/
    2057 void  Tria::MaxVy(double* pmaxvy, bool process_units){
    2058 
    2059         /*Get maximum:*/
    2060         double maxvy=this->inputs->Max(VyEnum);
    2061 
    2062         /*process units if requested: */
    2063         if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum);
    2064 
    2065         /*Assign output pointers:*/
    2066         *pmaxvy=maxvy;
    2067 
    2068 }
    2069 /*}}}*/
    2070 /*FUNCTION Tria::MaxVz{{{1*/
    2071 void  Tria::MaxVz(double* pmaxvz, bool process_units){
    2072 
    2073         /*Get maximum:*/
    2074         double maxvz=this->inputs->Max(VzEnum);
    2075 
    2076         /*process units if requested: */
    2077         if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum);
    2078 
    2079         /*Assign output pointers:*/
    2080         *pmaxvz=maxvz;
    2081 }
    2082 /*}}}*/
    20831928/*FUNCTION Tria::MigrateGroundingLine{{{1*/
    20841929void  Tria::MigrateGroundingLine(void){
    … …  
    21562001}
    21572002/*}}}*/
    2158 /*FUNCTION Tria::MinVel{{{1*/
    2159 void  Tria::MinVel(double* pminvel, bool process_units){
    2160 
    2161         /*Get minimum:*/
    2162         double minvel=this->inputs->Min(VelEnum);
    2163 
    2164         /*process units if requested: */
    2165         if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum);
    2166 
    2167         /*Assign output pointers:*/
    2168         *pminvel=minvel;
    2169 }
    2170 /*}}}*/
    2171 /*FUNCTION Tria::MinVx{{{1*/
    2172 void  Tria::MinVx(double* pminvx, bool process_units){
    2173 
    2174         /*Get minimum:*/
    2175         double minvx=this->inputs->Min(VxEnum);
    2176 
    2177         /*process units if requested: */
    2178         if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum);
    2179 
    2180         /*Assign output pointers:*/
    2181         *pminvx=minvx;
    2182 }
    2183 /*}}}*/
    2184 /*FUNCTION Tria::MinVy{{{1*/
    2185 void  Tria::MinVy(double* pminvy, bool process_units){
    2186 
    2187         /*Get minimum:*/
    2188         double minvy=this->inputs->Min(VyEnum);
    2189 
    2190         /*process units if requested: */
    2191         if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum);
    2192 
    2193         /*Assign output pointers:*/
    2194         *pminvy=minvy;
    2195 }
    2196 /*}}}*/
    2197 /*FUNCTION Tria::MinVz{{{1*/
    2198 void  Tria::MinVz(double* pminvz, bool process_units){
    2199 
    2200         /*Get minimum:*/
    2201         double minvz=this->inputs->Min(VzEnum);
    2202 
    2203         /*process units if requested: */
    2204         if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum);
    2205 
    2206         /*Assign output pointers:*/
    2207         *pminvz=minvz;
    2208 }
    2209 /*}}}*/
    22102003/*FUNCTION Tria::MyRank {{{1*/
    22112004int    Tria::MyRank(void){
    22122005        extern int my_rank;
    22132006        return my_rank;
    2214 }
    2215 /*}}}*/
    2216 /*FUNCTION Tria::NodalValue {{{1*/
    2217 int    Tria::NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units){
    2218 
    2219         int i;
    2220         int found=0;
    2221         double value;
    2222         Input* data=NULL;
    2223         GaussTria *gauss                            = NULL;
    2224 
    2225         /*First, serarch the input: */
    2226         data=inputs->GetInput(natureofdataenum);
    2227 
    2228         /*figure out if we have the vertex id: */
    2229         found=0;
    2230         for(i=0;i<NUMVERTICES;i++){
    2231                 if(index==nodes[i]->GetVertexId()){
    2232                         /*Do we have natureofdataenum in our inputs? :*/
    2233                         if(data){
    2234                                 /*ok, we are good. retrieve value of input at vertex :*/
    2235                                 gauss=new GaussTria(); gauss->GaussVertex(i);
    2236                                 data->GetParameterValue(&value,gauss);
    2237                                 found=1;
    2238                                 break;
    2239                         }
    2240                 }
    2241         }
    2242 
    2243         if(found)*pvalue=value;
    2244         return found;
    22452007}
    22462008/*}}}*/
    … …  
    23462108                elementresult->ProcessUnits(this->parameters);
    23472109        }
    2348 }
    2349 /*}}}*/
    2350 /*FUNCTION Tria::RheologyBbarx{{{1*/
    2351 double Tria::RheologyBbarx(void){
    2352 
    2353         return this->matice->GetBbar();
    2354 
    23552110}
    23562111/*}}}*/
    … …  
    25982353
    25992354        /*Get for Vx and Vy, the max of abs value: */
     2355        #ifdef _HAVE_RESPONSES_
    26002356        this->MaxAbsVx(&maxabsvx,false);
    26012357        this->MaxAbsVy(&maxabsvy,false);
     2358        #else
     2359                _error_("ISSM was not compiled with responses compiled in, exiting!");
     2360        #endif
    26022361
    26032362        /* Get node coordinates and dof list: */
    … …  
    28482607}
    28492608/*}}}*/
     2609
     2610#ifdef _HAVE_RESPONSES_
     2611/*FUNCTION Tria::MassFlux {{{1*/
     2612double Tria::MassFlux( double* segment,bool process_units){
     2613
     2614        const int    numdofs=2;
     2615
     2616        int        i;
     2617        double     mass_flux=0;
     2618        double     xyz_list[NUMVERTICES][3];
     2619        double     normal[2];
     2620        double     length,rho_ice;
     2621        double     x1,y1,x2,y2,h1,h2;
     2622        double     vx1,vx2,vy1,vy2;
     2623        GaussTria* gauss_1=NULL;
     2624        GaussTria* gauss_2=NULL;
     2625
     2626        /*Get material parameters :*/
     2627        rho_ice=matpar->GetRhoIce();
     2628
     2629        /*First off, check that this segment belongs to this element: */
     2630        if ((int)*(segment+4)!=this->id)_error_("%s%i%s%i","error message: segment with id ",(int)*(segment+4)," does not belong to element with id:",this->id);
     2631
     2632        /*Recover segment node locations: */
     2633        x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
     2634
     2635        /*Get xyz list: */
     2636        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     2637
     2638        /*get area coordinates of 0 and 1 locations: */
     2639        gauss_1=new GaussTria();
     2640        gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
     2641        gauss_2=new GaussTria();
     2642        gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
     2643
     2644        normal[0]=cos(atan2(x1-x2,y2-y1));
     2645        normal[1]=sin(atan2(x1-x2,y2-y1));
     2646
     2647        length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2));
     2648
     2649        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     2650        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     2651        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     2652
     2653        thickness_input->GetParameterValue(&h1, gauss_1);
     2654        thickness_input->GetParameterValue(&h2, gauss_2);
     2655        vx_input->GetParameterValue(&vx1,gauss_1);
     2656        vx_input->GetParameterValue(&vx2,gauss_2);
     2657        vy_input->GetParameterValue(&vy1,gauss_1);
     2658        vy_input->GetParameterValue(&vy2,gauss_2);
     2659
     2660        mass_flux= rho_ice*length*( 
     2661                                (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
     2662                                (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
     2663                                );
     2664
     2665        /*Process units: */
     2666        mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum);
     2667
     2668        /*clean up and return:*/
     2669        delete gauss_1;
     2670        delete gauss_2;
     2671        return mass_flux;
     2672}
     2673/*}}}*/
     2674/*FUNCTION Tria::MaxAbsVx{{{1*/
     2675void  Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){
     2676
     2677        /*Get maximum:*/
     2678        double maxabsvx=this->inputs->MaxAbs(VxEnum);
     2679
     2680        /*process units if requested: */
     2681        if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum);
     2682
     2683        /*Assign output pointers:*/
     2684        *pmaxabsvx=maxabsvx;
     2685}
     2686/*}}}*/
     2687/*FUNCTION Tria::MaxAbsVy{{{1*/
     2688void  Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){
     2689
     2690        /*Get maximum:*/
     2691        double maxabsvy=this->inputs->MaxAbs(VyEnum);
     2692
     2693        /*process units if requested: */
     2694        if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum);
     2695
     2696        /*Assign output pointers:*/
     2697        *pmaxabsvy=maxabsvy;
     2698}
     2699/*}}}*/
     2700/*FUNCTION Tria::MaxAbsVz{{{1*/
     2701void  Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){
     2702
     2703        /*Get maximum:*/
     2704        double maxabsvz=this->inputs->MaxAbs(VzEnum);
     2705
     2706        /*process units if requested: */
     2707        if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum);
     2708
     2709        /*Assign output pointers:*/
     2710        *pmaxabsvz=maxabsvz;
     2711}
     2712/*}}}*/
     2713/*FUNCTION Tria::MaxVel{{{1*/
     2714void  Tria::MaxVel(double* pmaxvel, bool process_units){
     2715
     2716        /*Get maximum:*/
     2717        double maxvel=this->inputs->Max(VelEnum);
     2718
     2719        /*process units if requested: */
     2720        if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum);
     2721
     2722        /*Assign output pointers:*/
     2723        *pmaxvel=maxvel;
     2724}
     2725/*}}}*/
     2726/*FUNCTION Tria::MaxVx{{{1*/
     2727void  Tria::MaxVx(double* pmaxvx, bool process_units){
     2728
     2729        /*Get maximum:*/
     2730        double maxvx=this->inputs->Max(VxEnum);
     2731
     2732        /*process units if requested: */
     2733        if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum);
     2734
     2735        /*Assign output pointers:*/
     2736        *pmaxvx=maxvx;
     2737}
     2738/*}}}*/
     2739/*FUNCTION Tria::MaxVy{{{1*/
     2740void  Tria::MaxVy(double* pmaxvy, bool process_units){
     2741
     2742        /*Get maximum:*/
     2743        double maxvy=this->inputs->Max(VyEnum);
     2744
     2745        /*process units if requested: */
     2746        if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum);
     2747
     2748        /*Assign output pointers:*/
     2749        *pmaxvy=maxvy;
     2750
     2751}
     2752/*}}}*/
     2753/*FUNCTION Tria::MaxVz{{{1*/
     2754void  Tria::MaxVz(double* pmaxvz, bool process_units){
     2755
     2756        /*Get maximum:*/
     2757        double maxvz=this->inputs->Max(VzEnum);
     2758
     2759        /*process units if requested: */
     2760        if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum);
     2761
     2762        /*Assign output pointers:*/
     2763        *pmaxvz=maxvz;
     2764}
     2765/*}}}*/
     2766/*FUNCTION Tria::MinVel{{{1*/
     2767void  Tria::MinVel(double* pminvel, bool process_units){
     2768
     2769        /*Get minimum:*/
     2770        double minvel=this->inputs->Min(VelEnum);
     2771
     2772        /*process units if requested: */
     2773        if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum);
     2774
     2775        /*Assign output pointers:*/
     2776        *pminvel=minvel;
     2777}
     2778/*}}}*/
     2779/*FUNCTION Tria::MinVx{{{1*/
     2780void  Tria::MinVx(double* pminvx, bool process_units){
     2781
     2782        /*Get minimum:*/
     2783        double minvx=this->inputs->Min(VxEnum);
     2784
     2785        /*process units if requested: */
     2786        if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum);
     2787
     2788        /*Assign output pointers:*/
     2789        *pminvx=minvx;
     2790}
     2791/*}}}*/
     2792/*FUNCTION Tria::MinVy{{{1*/
     2793void  Tria::MinVy(double* pminvy, bool process_units){
     2794
     2795        /*Get minimum:*/
     2796        double minvy=this->inputs->Min(VyEnum);
     2797
     2798        /*process units if requested: */
     2799        if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum);
     2800
     2801        /*Assign output pointers:*/
     2802        *pminvy=minvy;
     2803}
     2804/*}}}*/
     2805/*FUNCTION Tria::MinVz{{{1*/
     2806void  Tria::MinVz(double* pminvz, bool process_units){
     2807
     2808        /*Get minimum:*/
     2809        double minvz=this->inputs->Min(VzEnum);
     2810
     2811        /*process units if requested: */
     2812        if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum);
     2813
     2814        /*Assign output pointers:*/
     2815        *pminvz=minvz;
     2816}
     2817/*}}}*/
     2818/*FUNCTION Tria::NodalValue {{{1*/
     2819int    Tria::NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units){
     2820
     2821        int i;
     2822        int found=0;
     2823        double value;
     2824        Input* data=NULL;
     2825        GaussTria *gauss                            = NULL;
     2826
     2827        /*First, serarch the input: */
     2828        data=inputs->GetInput(natureofdataenum);
     2829
     2830        /*figure out if we have the vertex id: */
     2831        found=0;
     2832        for(i=0;i<NUMVERTICES;i++){
     2833                if(index==nodes[i]->GetVertexId()){
     2834                        /*Do we have natureofdataenum in our inputs? :*/
     2835                        if(data){
     2836                                /*ok, we are good. retrieve value of input at vertex :*/
     2837                                gauss=new GaussTria(); gauss->GaussVertex(i);
     2838                                data->GetParameterValue(&value,gauss);
     2839                                found=1;
     2840                                break;
     2841                        }
     2842                }
     2843        }
     2844
     2845        if(found)*pvalue=value;
     2846        return found;
     2847}
     2848/*}}}*/
     2849/*FUNCTION Tria::RheologyBbarx{{{1*/
     2850double Tria::RheologyBbarx(void){
     2851
     2852        return this->matice->GetBbar();
     2853
     2854}
     2855/*}}}*/
     2856#endif
    28502857
    28512858#ifdef _HAVE_DIAGNOSTIC_
Note: See TracChangeset for help on using the changeset viewer.