Changeset 15060 for issm/trunk-jpl/src


Ignore:
Timestamp:
05/21/13 15:11:26 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: moving BasisInt to generic functions + some cleanup

Location:
issm/trunk-jpl/src/c/classes/Elements
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r15015 r15060  
    4141                virtual void   CreatePVector(Vector<IssmDouble>* pf)=0;
    4242                virtual void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0;
     43                virtual void   BasisIntegral(Vector<IssmDouble>* basisg)=0;
    4344                virtual void   GetSolutionFromInputs(Vector<IssmDouble>* solution)=0;
    4445                virtual int    GetNodeIndex(Node* node)=0;
    … …  
    135136                #ifdef _HAVE_HYDROLOGY_
    136137                virtual void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode)=0;
    137                 virtual void BasisIntegral(Vector<IssmDouble>* basisg)=0;
    138138                virtual void GetHydrologyTransfer(Vector<IssmDouble>* transfer)=0;
    139139                #endif

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

    r15042 r15060  
    802802        /*return output*/
    803803        return penta;
     804}
     805/*}}}*/
     806/*FUNCTION Penta::BasisIntegral {{{*/
     807void Penta::BasisIntegral(Vector<IssmDouble>* basisg){
     808
     809        /*Constants*/
     810        const int    numdof=NDOF1*NUMVERTICES;
     811
     812        /*Intermediaries */
     813        IssmDouble Jdet;
     814        IssmDouble xyz_list[NUMVERTICES][3];
     815        IssmDouble basis[numdof];
     816        IssmDouble basisint[numdof]={0.};
     817        int       *doflist=NULL;
     818        GaussPenta* gauss=NULL;
     819
     820        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     821        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     822
     823        /* Start looping on the number of gaussian points: */
     824        gauss=new GaussPenta(2,2);
     825        for(int ig=gauss->begin();ig<gauss->end();ig++){
     826
     827                gauss->GaussPoint(ig);
     828
     829                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     830                GetNodalFunctionsP1(&basis[0], gauss);
     831
     832                for(int i=0;i<numdof;i++) basisint[i]+=Jdet*gauss->weight*basis[i];
     833        }
     834
     835        basisg->SetValues(numdof,doflist,&basisint[0],ADD_VAL);
     836
     837        /*Clean up and return*/
     838        delete gauss;
     839        xDelete<int>(doflist);
    804840}
    805841/*}}}*/
    … …  
    20732109void  Penta::InputUpdateFromVector(IssmDouble* vector, int name, int type){
    20742110
     2111        const int   numdof         = NDOF1 *NUMVERTICES;
     2112        int        *doflist        = NULL;
     2113        IssmDouble  values[numdof];
     2114
    20752115        /*Check that name is an element input*/
    20762116        if (!IsInput(name)) return;
    20772117
    2078         /*Penta update B in InputUpdateFromSolutionThermal, so don't look for B update here.*/
    2079 
    20802118        switch(type){
    2081 
    2082                 case VertexEnum:
    2083                         {
    2084 
    2085                                 /*New PentaVertexInpu*/
    2086                                 IssmDouble values[6];
    2087 
    2088                                 /*Get values on the 6 vertices*/
    2089                                 for (int i=0;i<6;i++){
    2090                                         values[i]=vector[this->vertices[i]->Pid()];
    2091                                 }
    2092 
    2093                                 /*update input*/
    2094                                 if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
    2095                                         material->inputs->AddInput(new PentaP1Input(name,values));
    2096                                 }
    2097                                 else{
    2098                                         this->inputs->AddInput(new PentaP1Input(name,values));
    2099                                 }
    2100                                 return;
    2101                                 break;
     2119                case VertexPIdEnum:
     2120                        for (int i=0;i<NUMVERTICES;i++){
     2121                                values[i]=vector[this->vertices[i]->Pid()];
    21022122                        }
     2123                        /*update input*/
     2124                        if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
     2125                                material->inputs->AddInput(new PentaP1Input(name,values));
     2126                        }
     2127                        else{
     2128                                this->inputs->AddInput(new PentaP1Input(name,values));
     2129                        }
     2130                        return;
     2131
     2132                case VertexSIdEnum:
     2133                        for (int i=0;i<NUMVERTICES;i++){
     2134                                values[i]=vector[this->vertices[i]->Sid()];
     2135                        }
     2136                        /*update input*/
     2137                        if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
     2138                                material->inputs->AddInput(new PentaP1Input(name,values));
     2139                        }
     2140                        else{
     2141                                this->inputs->AddInput(new PentaP1Input(name,values));
     2142                        }
     2143                        return;
     2144
     2145                case NodesEnum:
     2146                        /*Get dof list: */
     2147                        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     2148
     2149                        /*Use the dof list to index into the vector: */
     2150                        for(int i=0;i<numdof;i++){
     2151                                values[i]=vector[doflist[i]];
     2152                                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
     2153                        }
     2154                        /*Add input to the element: */
     2155                        this->inputs->AddInput(new PentaP1Input(name,values));
     2156
     2157                        /*Free ressources:*/
     2158                        xDelete<int>(doflist);
     2159                        return;
    21032160
    21042161                default:
    2105 
    21062162                        _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
    21072163        }
    … …  
    21482204                                name==InversionVyObsEnum ||
    21492205                                name==InversionVzObsEnum ||
     2206                                name==BasisIntegralEnum ||
    21502207                                name==TemperatureEnum ||
    21512208                                name==EnthalpyEnum ||
    … …  
    45144571        const int    numdof=NDOF1*NUMVERTICES;
    45154572
    4516         bool   converged;
    4517         int    i,rheology_law;
    4518         IssmDouble xyz_list[NUMVERTICES][3];
    4519         IssmDouble values[numdof];
    4520         IssmDouble B[numdof];
    4521         IssmDouble B_average,s_average;
    4522         int*   doflist=NULL;
    4523         //IssmDouble pressure[numdof];
     4573        bool        converged;
     4574        int         i,rheology_law;
     4575        IssmDouble  xyz_list[NUMVERTICES][3];
     4576        IssmDouble  values[numdof];
     4577        IssmDouble  B[numdof];
     4578        IssmDouble  B_average,s_average;
     4579        int        *doflist = NULL;
     4580        IssmDouble pressure[numdof];
    45244581
    45254582        /*Get dof list: */
    … …  
    93189375}
    93199376/*}}}*/
    9320 /*FUNCTION Penta::BasisIntegral {{{*/
    9321 void Penta::BasisIntegral(Vector<IssmDouble>* basisg){
    9322         _error_("Hydrological stuff not suported in Penta");
    9323 }
    9324 /*}}}*/
    9325 
    9326 /*}}}*/
    93279377/*FUNCTION Penta::GetHydrologyTransfer{{{*/
    93289378void  Penta::GetHydrologyTransfer(Vector<IssmDouble>* transfer){
    … …  
    93309380}
    93319381/*}}}*/
    9332 
    93339382#endif

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

    r15049 r15060  
    265265                void           InputUpdateFromSolutionDiagnosticVert( IssmDouble* solutiong);
    266266                void           InputUpdateFromSolutionDiagnosticStokes( IssmDouble* solutiong);
     267                void           BasisIntegral(Vector<IssmDouble>* gbasis);
    267268                void             GetSolutionFromInputsDiagnosticHoriz(Vector<IssmDouble>* solutiong);
    268269                void             GetSolutionFromInputsDiagnosticHutter(Vector<IssmDouble>* solutiong);
    … …  
    307308                void    CreateHydrologyWaterVelocityInput(void);
    308309                void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    309                 void    BasisIntegral(Vector<IssmDouble>* gbasis);
    310310                void    GetHydrologyTransfer(Vector<IssmDouble>* transfer);
    311311                #endif

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

    r15054 r15060  
    10991099        _assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0);
    11001100        return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2;
     1101}
     1102/*}}}*/
     1103/*FUNCTION Tria::BasisIntegral {{{*/
     1104void Tria::BasisIntegral(Vector<IssmDouble>* basisg){
     1105
     1106        /*Constants*/
     1107        const int    numdof=NDOF1*NUMVERTICES;
     1108
     1109        /*Intermediaries */
     1110        IssmDouble Jdet;
     1111        IssmDouble xyz_list[NUMVERTICES][3];
     1112        IssmDouble basis[numdof];
     1113        IssmDouble basisint[numdof]={0.};
     1114        int       *doflist=NULL;
     1115        GaussTria* gauss=NULL;
     1116
     1117        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     1118        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     1119
     1120        /* Start looping on the number of gaussian points: */
     1121        gauss=new GaussTria(2);
     1122        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1123
     1124                gauss->GaussPoint(ig);
     1125
     1126                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     1127                GetNodalFunctions(&basis[0], gauss);
     1128
     1129                for(int i=0;i<numdof;i++) basisint[i]+=Jdet*gauss->weight*basis[i];
     1130        }
     1131
     1132        basisg->SetValues(numdof,doflist,&basisint[0],ADD_VAL);
     1133
     1134        /*Clean up and return*/
     1135        delete gauss;
     1136        xDelete<int>(doflist);
    11011137}
    11021138/*}}}*/
    … …  
    19381974        switch(type){
    19391975        case VertexPIdEnum:
    1940                 /*Get values on the 3 vertices*/
    1941                 for (int i=0;i<3;i++){
     1976                for (int i=0;i<NUMVERTICES;i++){
    19421977                        values[i]=vector[this->vertices[i]->Pid()];
    19431978                }
    … …  
    19521987
    19531988        case VertexSIdEnum:
    1954                 /*Get values on the 3 vertices*/
    1955                 for (int i=0;i<3;i++){
     1989                for (int i=0;i<NUMVERTICES;i++){
    19561990                        values[i]=vector[this->vertices[i]->Sid()];
    19571991                }
    … …  
    64876521}
    64886522/*}}}*/
    6489 /*FUNCTION Tria::BasisIntegral {{{*/
    6490 void Tria::BasisIntegral(Vector<IssmDouble>* basisg){
    6491 
    6492         /*Constants*/
    6493         const int    numdof=NDOF1*NUMVERTICES;
    6494 
    6495         /*Intermediaries */
    6496         IssmDouble Jdet;
    6497         IssmDouble xyz_list[NUMVERTICES][3];
    6498         IssmDouble basis[numdof];
    6499         IssmDouble basisint[numdof]={0.};
    6500         int       *doflist=NULL;
    6501         GaussTria* gauss=NULL;
    6502 
    6503         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    6504         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6505 
    6506         /* Start looping on the number of gaussian points: */
    6507         gauss=new GaussTria(2);
    6508         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6509 
    6510                 gauss->GaussPoint(ig);
    6511 
    6512                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6513                 GetNodalFunctions(&basis[0], gauss);
    6514 
    6515                 for(int i=0;i<numdof;i++) basisint[i]+=Jdet*gauss->weight*basis[i];
    6516         }
    6517 
    6518         basisg->SetValues(numdof,doflist,&basisint[0],ADD_VAL);
    6519 
    6520         /*Clean up and return*/
    6521         delete gauss;
    6522         xDelete<int>(doflist);
    6523 }
    6524 /*}}}*/
    65256523/*FUNCTION Tria::GetHydrologyTransfer{{{*/
    65266524void  Tria::GetHydrologyTransfer(Vector<IssmDouble>* transfer){
    … …  
    65826580
    65836581/*}}}*/
    6584 
    65856582#endif
    65866583

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

    r15054 r15060  
    198198                ElementVector* CreatePVectorSlope(void);
    199199                IssmDouble     GetArea(void);
     200                void           BasisIntegral(Vector<IssmDouble>* gbasis);
    200201                int            GetElementType(void);
    201202                void             GetDofList(int** pdoflist,int approximation_enum,int setenum);
    … …  
    255256                void      InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution);
    256257                void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    257                 void    BasisIntegral(Vector<IssmDouble>* gbasis);
    258258                void    GetHydrologyTransfer(Vector<IssmDouble>* transfer);
    259259                #endif
Note: See TracChangeset for help on using the changeset viewer.