Changeset 5633


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

more GaussTria

Location:
issm/trunk/src/c/objects
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r5578 r5633  
    3535                virtual bool   GetShelf()=0;
    3636                virtual bool   GetOnBed()=0;
    37                 virtual void   GetThicknessList(double* thickness_list)=0;
     37                virtual void   GetParameterListOnVertices(double* pvalue,int enumtype)=0;
     38                virtual void   GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue)=0;
    3839                virtual void   GetParameterValue(double* pvalue,Node* node,int enumtype)=0;
    3940                virtual void   GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in)=0;
    40                 virtual void   GetBedList(double* bed_list)=0;
    4141                virtual void   Gradj(Vec gradient,int control_type)=0;
    4242                virtual void   GradjDrag(Vec gradient)=0;

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

    r5624 r5633  
    943943}
    944944/*}}}*/
    945 /*FUNCTION Penta::GetBedList {{{1*/
    946 void Penta::GetBedList(double* bedlist){
    947 
    948         const int numgrids=6;
    949         double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
    950        
    951         inputs->GetParameterValues(bedlist,&gaussgrids[0][0],6,BedEnum);
    952 
    953 }
    954 /*}}}*/
    955945/*FUNCTION Penta::GetMatPar {{{1*/
    956946void* Penta::GetMatPar(){
    … …  
    10241014                ISSMERROR("analysis: %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
    10251015        }
    1026 }
    1027 /*}}}*/
    1028 /*FUNCTION Penta::GetThicknessList {{{1*/
    1029 void Penta::GetThicknessList(double* thickness_list){
    1030 
    1031         const int numgrids=6;
    1032         double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
    1033         inputs->GetParameterValues(thickness_list,&gaussgrids[0][0],6,ThicknessEnum);
    1034 
    10351016}
    10361017/*}}}*/
    … …  
    48354816}
    48364817/*}}}*/
     4818/*FUNCTION Penta::GetParameterListOnVertices(double* pvalue,int enumtype) {{{1*/
     4819void Penta::GetParameterListOnVertices(double* pvalue,int enumtype){
     4820
     4821        /*Intermediaries*/
     4822        const int  numvertices = 6;
     4823        double     value[numvertices];
     4824        double     gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
     4825
     4826        /*Recover input*/
     4827        Input* input=inputs->GetInput(enumtype);
     4828        if (!input) ISSMERROR("Input %s not found in element",EnumToString(enumtype));
     4829
     4830        /*Checks in debugging mode*/
     4831        ISSMASSERT(pvalue);
     4832
     4833        /* Start looping on the number of vertices: */
     4834        for (int iv=0;iv<numvertices;iv++){
     4835                input->GetParameterValue(&pvalue[iv],&gauss[iv][0]);
     4836        }
     4837
     4838        /*clean-up*/
     4839}
     4840/*}}}*/
     4841/*FUNCTION Penta::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue) {{{1*/
     4842void Penta::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue){
     4843
     4844        /*Intermediaries*/
     4845        const int  numvertices = 6;
     4846        double     value[numvertices];
     4847        double     gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
     4848
     4849        /*Recover input*/
     4850        Input* input=inputs->GetInput(enumtype);
     4851
     4852        /*Checks in debugging mode*/
     4853        ISSMASSERT(pvalue);
     4854
     4855        /* Start looping on the number of vertices: */
     4856        if (input)
     4857         for (int iv=0;iv<numvertices;iv++) input->GetParameterValue(&pvalue[iv],&gauss[iv][0]);
     4858        else
     4859         for (int iv=0;iv<numvertices;iv++) pvalue[iv]=defaultvalue;
     4860
     4861        /*clean-up*/
     4862        delete gauss;
     4863}
     4864/*}}}*/
    48374865/*FUNCTION Penta::GetParameterValue(double* pvalue,Node* node,int enumtype) {{{1*/
    48384866void Penta::GetParameterValue(double* pvalue,Node* node,int enumtype){

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

    r5590 r5633  
    7979                void   CreatePVector(Vec pg);
    8080                void   DeleteResults(void);
    81                 void   GetBedList(double* bed_list);
    8281                void*  GetMatPar();
    8382                void   GetNodes(void** nodes);
    … …  
    8584                bool   GetShelf();
    8685                void   GetSolutionFromInputs(Vec solution);
    87                 void   GetThicknessList(double* thickness_list);
    8886                void   GetVectorFromInputs(Vec vector,int NameEnum);
    8987                void   Gradj(Vec gradient,int control_type);
    … …  
    160158                void      GetNodalFunctionsDerivativesStokes(double* dh1dh7,double* xyz_list, double* gauss_coord);
    161159                void      GetNodalFunctionsStokes(double* l1l7, double* gauss_coord);
     160                void    GetParameterListOnVertices(double* pvalue,int enumtype);
     161                void    GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue);
    162162                void      GetParameterValue(double* pvalue,Node* node,int enumtype);
    163163                void      GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in);

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

    r5631 r5633  
    450450       
    451451        const int numvertices=3;
    452         double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    453452        int i,j;
    454453
    … …  
    477476                                       
    478477                                        /*retrieve inputs: */
    479                                         inputs->GetParameterValues(&bed[0],&gaussgrids[0][0],3,BedEnum);
    480                                         inputs->GetParameterValues(&surface[0],&gaussgrids[0][0],3,SurfaceEnum);
     478                                        GetParameterListOnVertices(&bed[0],BedEnum);
     479                                        GetParameterListOnVertices(&surface[0],SurfaceEnum);
    481480
    482481                                        /*build new thickness: */
    … …  
    594593        double thickness[numvertices];
    595594        double rho_ice,g;
    596         double gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    597595       
    598596        /*Get dof list on which we will plug the pressure values: */
    … …  
    602600        rho_ice=matpar->GetRhoIce();
    603601        g=matpar->GetG();
    604 
    605         /*recover value of thickness at gauss points (0,0,1), (0,1,0),(1,0,0): */
    606         inputs->GetParameterValues(&thickness[0],&gauss[0][0],3,ThicknessEnum);
    607 
    608         for(i=0;i<numvertices;i++){
    609                 pressure[i]=rho_ice*g*thickness[i];
    610         }
     602        GetParameterListOnVertices(&thickness[0],ThicknessEnum);
     603
     604        for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
    611605
    612606        /*plug local pressure values into global pressure vector: */
    … …  
    849843}
    850844/*}}}*/
    851 /*FUNCTION Tria::GetBedList {{{1*/
    852 void  Tria::GetBedList(double* bedlist){
    853 
    854         Input     *bed_input = NULL;
    855         GaussTria *gauss     = NULL;
    856 
    857         /*Recover input*/
    858         bed_input=inputs->GetInput(BedEnum);
    859 
    860         /*Checks in debugging mode*/
    861         ISSMASSERT(bedlist);
    862         ISSMASSERT(bed_input);
    863 
    864         /* Start looping on the number of vertices: */
    865         gauss=new GaussTria();
    866         for (int iv=0;iv<3;iv++){
    867                 gauss->GaussVertex(iv);
    868                 bed_input->GetParameterValue(&bedlist[iv],gauss);
    869         }
    870 
    871         /*clean-up*/
    872         delete gauss;
    873 
    874 }
    875 /*}}}*/
    876845/*FUNCTION Tria::GetMatPar {{{1*/
    877846void* Tria::GetMatPar(){
    … …  
    923892        if (analysis_type==DiagnosticHorizAnalysisEnum)
    924893         GetSolutionFromInputsDiagnosticHoriz(solution);
    925         else if (analysis_type==AdjointHorizAnalysisEnum)
    926          GetSolutionFromInputsAdjointHoriz(solution);
    927894        else if (analysis_type==DiagnosticHutterAnalysisEnum)
    928895         GetSolutionFromInputsDiagnosticHutter(solution);
    929896        else
    930897         ISSMERROR("analysis: %s not supported yet",EnumToString(analysis_type));
    931 
    932 }
    933 /*}}}*/
    934 /*FUNCTION Tria::GetThicknessList {{{1*/
    935 void Tria::GetThicknessList(double* thicknesslist){
    936 
    937         Input     *thickness_input = NULL;
    938         GaussTria *gauss     = NULL;
    939 
    940         /*Recover input*/
    941         thickness_input=inputs->GetInput(ThicknessEnum);
    942 
    943         /*Checks in debugging mode*/
    944         ISSMASSERT(thicknesslist);
    945         ISSMASSERT(thickness_input);
    946 
    947         /* Start looping on the number of vertices: */
    948         gauss=new GaussTria();
    949         for (int iv=0;iv<3;iv++){
    950                 gauss->GaussVertex(iv);
    951                 thickness_input->GetParameterValue(&thicknesslist[iv],gauss);
    952         }
    953 
    954         /*clean-up*/
    955         delete gauss;
    956898
    957899}
    … …  
    63386280}
    63396281/*}}}*/
     6282/*FUNCTION Tria::GetParameterListOnVertices(double* pvalue,int enumtype) {{{1*/
     6283void Tria::GetParameterListOnVertices(double* pvalue,int enumtype){
     6284
     6285        /*Intermediaries*/
     6286        const int  numvertices        = 3;
     6287        double     value[numvertices];
     6288        GaussTria *gauss              = NULL;
     6289
     6290        /*Recover input*/
     6291        Input* input=inputs->GetInput(enumtype);
     6292        if (!input) ISSMERROR("Input %s not found in element",EnumToString(enumtype));
     6293
     6294        /*Checks in debugging mode*/
     6295        ISSMASSERT(pvalue);
     6296
     6297        /* Start looping on the number of vertices: */
     6298        gauss=new GaussTria();
     6299        for (int iv=0;iv<3;iv++){
     6300                gauss->GaussVertex(iv);
     6301                input->GetParameterValue(&pvalue[iv],gauss);
     6302        }
     6303
     6304        /*clean-up*/
     6305        delete gauss;
     6306}
     6307/*}}}*/
     6308/*FUNCTION Tria::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue) {{{1*/
     6309void Tria::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue){
     6310
     6311        /*Intermediaries*/
     6312        const int  numvertices        = 3;
     6313        double     value[numvertices];
     6314        GaussTria *gauss              = NULL;
     6315
     6316        /*Recover input*/
     6317        Input* input=inputs->GetInput(enumtype);
     6318
     6319        /*Checks in debugging mode*/
     6320        ISSMASSERT(pvalue);
     6321
     6322        /* Start looping on the number of vertices: */
     6323        if (input){
     6324                gauss=new GaussTria();
     6325                for (int iv=0;iv<numvertices;iv++){
     6326                        gauss->GaussVertex(iv);
     6327                        input->GetParameterValue(&pvalue[iv],gauss);
     6328                }
     6329        }
     6330        else{
     6331                for (int iv=0;iv<numvertices;iv++) pvalue[iv]=defaultvalue;
     6332        }
     6333
     6334        /*clean-up*/
     6335        delete gauss;
     6336}
     6337/*}}}*/
    63406338/*FUNCTION Tria::GetParameterValue(double* pvalue,Node* node,int enumtype) {{{1*/
    63416339void Tria::GetParameterValue(double* pvalue,Node* node,int enumtype){
    … …  
    64956493
    64966494        int i;
    6497 
    64986495        const int    numvertices=3;
    64996496        const int    numdofpervertex=2;
    65006497        const int    numdof=numdofpervertex*numvertices;
    6501         double       gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    65026498        int*         doflist=NULL;
    65036499        double       values[numdof];
    65046500        double       vx;
    65056501        double       vy;
     6502        GaussTria*   gauss=NULL;
    65066503
    65076504        /*Get dof list: */
    65086505        GetDofList(&doflist);
    65096506
     6507        /*Get inputs*/
     6508        Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
     6509        Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
     6510
    65106511        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    65116512        /*P1 element only for now*/
     6513        gauss=new GaussTria();
    65126514        for(i=0;i<numvertices;i++){
    65136515
     6516                gauss->GaussVertex(i);
     6517
    65146518                /*Recover vx and vy*/
    6515                 inputs->GetParameterValue(&vx,&gauss[i][0],VxEnum);
    6516                 inputs->GetParameterValue(&vy,&gauss[i][0],VyEnum);
     6519                vx_input->GetParameterValue(&vx,gauss);
     6520                vy_input->GetParameterValue(&vy,gauss);
    65176521                values[i*numdofpervertex+0]=vx;
    65186522                values[i*numdofpervertex+1]=vy;
    … …  
    65236527
    65246528        /*Free ressources:*/
     6529        delete gauss;
    65256530        xfree((void**)&doflist);
    65266531
    65276532}
    65286533/*}}}*/
    6529 /*FUNCTION Tria::GetSolutionFromInputsAdjointHoriz{{{1*/
    6530 void  Tria::GetSolutionFromInputsAdjointHoriz(Vec solution){
     6534/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHutter{{{1*/
     6535void  Tria::GetSolutionFromInputsDiagnosticHutter(Vec solution){
    65316536
    65326537        int i;
    6533 
    65346538        const int    numvertices=3;
    65356539        const int    numdofpervertex=2;
    65366540        const int    numdof=numdofpervertex*numvertices;
    6537         double       gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    6538         int*         doflist=NULL;
    6539         double       values[numdof];
    6540         double       lambdax;
    6541         double       lambday;
    6542 
    6543         /*Get dof list: */
    6544         GetDofList(&doflist);
    6545 
    6546         /*Ok, we have lambdax and lambday in values, fill in lambdax and lambday arrays: */
    6547         /*P1 element only for now*/
    6548         for(i=0;i<numvertices;i++){
    6549 
    6550                 /*Recover lambdax and lambday*/
    6551                 inputs->GetParameterValue(&lambdax,&gauss[i][0],VxEnum);
    6552                 inputs->GetParameterValue(&lambday,&gauss[i][0],VyEnum);
    6553                 values[i*numdofpervertex+0]=lambdax;
    6554                 values[i*numdofpervertex+1]=lambday;
    6555         }
    6556 
    6557         /*Add value to global vector*/
    6558         VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
    6559 
    6560         /*Free ressources:*/
    6561         xfree((void**)&doflist);
    6562 
    6563 }
    6564 /*}}}*/
    6565 /*FUNCTION Tria::GetSolutionFromInputsDiagnosticHutter{{{1*/
    6566 void  Tria::GetSolutionFromInputsDiagnosticHutter(Vec solution){
    6567 
    6568         int i;
    6569 
    6570         const int    numvertices=3;
    6571         const int    numdofpervertex=2;
    6572         const int    numdof=numdofpervertex*numvertices;
    6573         double       gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    65746541        int*         doflist=NULL;
    65756542        double       values[numdof];
    … …  
    65776544        double       vy;
    65786545        int          dummy;
     6546        GaussTria*   gauss=NULL;
    65796547
    65806548        /*Get dof list: */
    65816549        GetDofList(&doflist);
    65826550
     6551        /*Get inputs*/
     6552        Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
     6553        Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
     6554
    65836555        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    65846556        /*P1 element only for now*/
     6557        gauss=new GaussTria();
    65856558        for(i=0;i<numvertices;i++){
    65866559
     6560                gauss->GaussVertex(i);
     6561
    65876562                /*Recover vx and vy*/
    6588                 inputs->GetParameterValue(&vx,&gauss[i][0],VxEnum);
    6589                 inputs->GetParameterValue(&vy,&gauss[i][0],VyEnum);
     6563                vx_input->GetParameterValue(&vx,gauss);
     6564                vy_input->GetParameterValue(&vy,gauss);
    65906565                values[i*numdofpervertex+0]=vx;
    65916566                values[i*numdofpervertex+1]=vy;
    … …  
    65966571
    65976572        /*Free ressources:*/
     6573        delete gauss;
    65986574        xfree((void**)&doflist);
    65996575
    … …  
    68926868        double       thickness[numvertices];
    68936869        double       rho_ice,g;
    6894         double       gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
    6895         Input*       vz_input=NULL;
    6896         double*      vz_ptr=NULL;
    68976870        int          dummy;
    68986871       
    … …  
    69126885
    69136886        /*Get Vz*/
    6914         vz_input=inputs->GetInput(VzEnum);
    6915         if (vz_input){
    6916                 if (vz_input->Enum()!=TriaVertexInputEnum){
    6917                         ISSMERROR("Cannot compute Vel as Vz is of type %s",EnumToString(vz_input->Enum()));
    6918                 }
    6919                 vz_input->GetValuesPtr(&vz_ptr,&dummy);
    6920                 for(i=0;i<numvertices;i++) vz[i]=vz_ptr[i];
    6921         }
    6922         else{
    6923                 for(i=0;i<numvertices;i++) vz[i]=0.0;
    6924         }
     6887        GetParameterListOnVertices(&vz[0],VzEnum,0);
    69256888
    69266889        /*Now Compute vel*/
    … …  
    69316894        rho_ice=matpar->GetRhoIce();
    69326895        g=matpar->GetG();
    6933         inputs->GetParameterValues(&thickness[0],&gauss[0][0],3,ThicknessEnum);
    6934        
    6935         for(i=0;i<numvertices;i++){
    6936                 pressure[i]=rho_ice*g*thickness[i];
    6937         }
     6896        GetParameterListOnVertices(&thickness[0],ThicknessEnum);
     6897        for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
    69386898
    69396899        /*Now, we have to move the previous Vx and Vy inputs  to old

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

    r5631 r5633  
    7575                void   CreateKMatrix(Mat Kgg);
    7676                void   CreatePVector(Vec pg);
    77                 void   GetBedList(double* bed_list);
    7877                void*  GetMatPar();
    7978                void   GetNodes(void** nodes);
    … …  
    8180                bool   GetShelf();
    8281                void   GetSolutionFromInputs(Vec solution);
    83                 void   GetThicknessList(double* thickness_list);
    8482                void   GetVectorFromInputs(Vec vector,int NameEnum);
    8583                void   Gradj(Vec gradient,int control_type);
    … …  
    155153                void      GetDofList(int** pdoflist,int approximation_enum=0);
    156154                void      GetDofList1(int* doflist);
     155                void    GetParameterListOnVertices(double* pvalue,int enumtype);
     156                void    GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue);
    157157                void      GetParameterValue(double* pvalue,Node* node,int enumtype);
    158158                void      GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype);
    159159                void      GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in);
    160160                void      GetSolutionFromInputsDiagnosticHoriz(Vec solution);
    161                 void      GetSolutionFromInputsAdjointHoriz(Vec solution);
    162161                void      GetSolutionFromInputsDiagnosticHutter(Vec solution);
    163162                void    GetStrainRate2d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input);

  • issm/trunk/src/c/objects/Loads/Icefront.cpp

    r5624 r5633  
    592592        }
    593593
    594         element->GetThicknessList(&thickness_list[0]);
    595         element->GetBedList(&bed_list[0]);
     594        element->GetParameterListOnVertices(&thickness_list[0],ThicknessEnum);
     595        element->GetParameterListOnVertices(&bed_list[0],BedEnum);
    596596
    597597        thickness_list_quad[0]=thickness_list[grid1];
    … …  
    731731        }
    732732
    733         element->GetThicknessList(&thickness_list[0]);
    734         element->GetBedList(&bed_list[0]);
     733        element->GetParameterListOnVertices(&thickness_list[0],ThicknessEnum);
     734        element->GetParameterListOnVertices(&bed_list[0],BedEnum);
    735735
    736736        thickness_list_quad[0]=thickness_list[grid1];
Note: See TracChangeset for help on using the changeset viewer.