Changeset 15714


Ignore:
Timestamp:
08/05/13 16:29:35 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: more changes so that we can have multiple finite elements

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

Legend:

Unmodified
Added
Removed

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

    r15712 r15714  
    48034803        const int    numdof=NDOF1*NUMVERTICES;
    48044804
    4805         int          i;
    4806         int*         doflist=NULL;
    4807         IssmDouble       values[numdof];
    4808         IssmDouble       enthalpy;
    4809         GaussPenta   *gauss=NULL;
     4805        int*        doflist=NULL;
     4806        IssmDouble  values[numdof];
     4807        IssmDouble  enthalpy;
     4808        GaussPenta *gauss=NULL;
    48104809
    48114810        /*Get dof list: */
    … …  
    48144813
    48154814        gauss=new GaussPenta();
    4816         for(i=0;i<NUMVERTICES;i++){
     4815        for(int i=0;i<NUMVERTICES;i++){
    48174816                /*Recover temperature*/
    48184817                gauss->GaussVertex(i);
    … …  
    53205319        int vnumnodes = this->NumberofNodesVelocity();
    53215320        int pnumnodes = this->NumberofNodesPressure();
    5322         int vnumdof   = vnumnodes*NDOF3;
    53235321
    53245322        /*Prepare coordinate system list*/
    … …  
    53295327        /*Initialize Element matrix and vectors*/
    53305328        ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
    5331         IssmDouble*    vbasis = xNew<IssmDouble>(vnumdof);
     5329        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    53325330
    53335331        /*Retrieve all inputs and parameters*/
    … …  
    77957793ElementMatrix* Penta::CreateKMatrixDiagnosticVertVolume(void){
    77967794
    7797         /*Constants*/
    7798         const int    numdof=NDOF1*NUMVERTICES;
    7799 
    78007795        /*Intermediaries */
    78017796        IssmDouble  Jdet;
    … …  
    78987893/*FUNCTION Penta::CreatePVectorCouplingSSAFSViscous {{{*/
    78997894ElementVector* Penta::CreatePVectorCouplingSSAFSViscous(void){
    7900 
    7901         /*Constants*/
    7902         const int   numdof=NUMVERTICES*NDOF4;
    79037895
    79047896        /*Intermediaries */
    … …  
    79607952/*FUNCTION Penta::CreatePVectorCouplingSSAFSFriction{{{*/
    79617953ElementVector* Penta::CreatePVectorCouplingSSAFSFriction(void){
    7962 
    7963         /*Constants*/
    7964         const int numdof=NUMVERTICES*NDOF4;
    79657954
    79667955        /*Intermediaries*/
    … …  
    80508039ElementVector* Penta::CreatePVectorCouplingHOFSViscous(void){
    80518040
    8052         /*Constants*/
    8053         const int   numdof=NUMVERTICES*NDOF4;
    8054 
    80558041        /*Intermediaries */
    80568042        int         i;
    … …  
    81118097/*FUNCTION Penta::CreatePVectorCouplingHOFSFriction{{{*/
    81128098ElementVector* Penta::CreatePVectorCouplingHOFSFriction(void){
    8113 
    8114         /*Constants*/
    8115         const int numdof=NUMVERTICES*NDOF4;
    81168099
    81178100        /*Intermediaries*/
    … …  
    82798262        /*Fetch number of nodes and dof for this finite element*/
    82808263        int numnodes = this->NumberofNodes(); _assert_(numnodes==6);
    8281         int numdof   = numnodes*NDOF2;
    82828264
    82838265        /*Initialize Element vector*/
    … …  
    84068388        /*Fetch number of nodes and dof for this finite element*/
    84078389        int numnodes = this->NumberofNodes();
    8408         int numdof   = numnodes*NDOF2;
    84098390
    84108391        /*Initialize vector*/
    … …  
    84798460        /*Fetch number of nodes and dof for this finite element*/
    84808461        int numnodes = this->NumberofNodes();
    8481         int numdof   = numnodes*NDOF2;
    84828462
    84838463        /*Initialize Element vector and other vectors*/
    … …  
    85868566        int vnumnodes = this->NumberofNodesVelocity();
    85878567        int pnumnodes = this->NumberofNodesPressure();
    8588         int vnumdof   = vnumnodes*NDOF3;
    85898568
    85908569        /*Prepare coordinate system list*/
    … …  
    85958574        /*Initialize Element matrix and vectors*/
    85968575        ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
    8597         IssmDouble*    vbasis = xNew<IssmDouble>(vnumdof);
     8576        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    85988577
    85998578        /*Retrieve all inputs and parameters*/
    … …  
    86428621/*FUNCTION Penta::PVectorGLSstabilization{{{*/
    86438622void Penta::PVectorGLSstabilization(ElementVector* pe){
    8644 
    8645         /*Constants*/
    8646         const int numdof=NDOF4*NUMVERTICES;
    86478623
    86488624        /*Intermediaries*/
    … …  
    87538729        int vnumnodes = this->NumberofNodesVelocity();
    87548730        int pnumnodes = this->NumberofNodesPressure();
    8755         int vnumdof   = vnumnodes*NDOF3;
    87568731
    87578732        /*Prepare coordinate system list*/
    … …  
    87628737        /*Initialize Element matrix and vectors*/
    87638738        ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
    8764         IssmDouble*    vbasis = xNew<IssmDouble>(vnumdof);
     8739        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    87658740
    87668741        /*Retrieve all inputs and parameters*/
    … …  
    88268801        int vnumnodes = this->NumberofNodesVelocity();
    88278802        int pnumnodes = this->NumberofNodesPressure();
    8828         int vnumdof   = vnumnodes*NDOF3;
    88298803
    88308804        /*Prepare coordinate system list*/
    … …  
    88358809        /*Initialize Element matrix and vectors*/
    88368810        ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
    8837         IssmDouble*    vbasis = xNew<IssmDouble>(vnumdof);
     8811        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    88388812
    88398813        /*Retrieve all inputs and parameters*/

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

    r15713 r15714  
    50045004ElementVector* Tria::CreatePVectorAdjointHoriz(void){
    50055005
    5006         /*Constants*/
    5007         const int    numdof=NDOF2*NUMVERTICES;
    5008 
    50095006        /*Intermediaries */
    5010         int         i,resp;
    5011         int        *responses=NULL;
    5012         int         num_responses;
     5007        int        i,resp;
     5008        int       *responses=NULL;
     5009        int        num_responses;
    50135010        IssmDouble Jdet;
    50145011        IssmDouble obs_velocity_mag,velocity_mag;
    … …  
    50165013        IssmDouble epsvel=2.220446049250313e-16;
    50175014        IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
    5018         IssmDouble scalex=0,scaley=0,scale=0,S=0;
     5015        IssmDouble scalex=0.,scaley=0.,scale=0.,S=0.;
    50195016        IssmDouble vx,vy,vxobs,vyobs,weight;
    50205017        IssmDouble xyz_list[NUMVERTICES][3];
    5021         IssmDouble basis[3];
    5022         GaussTria* gauss=NULL;
     5018
     5019        /*Fetch number of nodes and dof for this finite element*/
     5020        int numnodes = this->NumberofNodes();
     5021        int numdof   = numnodes*NDOF2;
    50235022
    50245023        /*Initialize Element vector*/
    5025         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     5024        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
     5025        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    50265026
    50275027        /*Retrieve all inputs and parameters*/
    … …  
    50435043
    50445044        /* Start  looping on the number of gaussian points: */
    5045         gauss=new GaussTria(4);
     5045        GaussTria* gauss=new GaussTria(4);
    50465046        for(int ig=gauss->begin();ig<gauss->end();ig++){
    50475047
    … …  
    50745074                                         *        du     obs
    50755075                                         */
    5076                                         for (i=0;i<NUMVERTICES;i++){
     5076                                        for(i=0;i<numnodes;i++){
    50775077                                                dux=vxobs-vx;
    50785078                                                duy=vyobs-vy;
    … …  
    50935093                                         *               obs
    50945094                                         */
    5095                                         for (i=0;i<NUMVERTICES;i++){
     5095                                        for(i=0;i<numnodes;i++){
    50965096                                                scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
    50975097                                                scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
    … …  
    51145114                                         *           
    51155115                                         */
    5116                                         for (i=0;i<NUMVERTICES;i++){
     5116                                        for(i=0;i<numnodes;i++){
    51175117                                                velocity_mag    =sqrt(pow(vx,   2)+pow(vy,   2))+epsvel;
    51185118                                                obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
    … …  
    51345134                                         *        du      S  2 sqrt(...)           obs
    51355135                                         */
    5136                                         for (i=0;i<NUMVERTICES;i++){
     5136                                        for(i=0;i<numnodes;i++){
    51375137                                                scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
    51385138                                                dux=scale*(vxobs-vx);
    … …  
    51525152                                         *        du                         |u| + eps  |u|                           u + eps
    51535153                                         */
    5154                                         for (i=0;i<NUMVERTICES;i++){
     5154                                        for(i=0;i<numnodes;i++){
    51555155                                                dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
    51565156                                                duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
    … …  
    51815181
    51825182        /*Clean up and return*/
     5183        xDelete<IssmDouble>(basis);
     5184        xDelete<int>(responses);
    51835185        delete gauss;
    5184         xDelete<int>(responses);
    51855186        return pe;
    51865187}
    … …  
    64946495        IssmDouble xyz_list[NUMVERTICES][3];
    64956496
    6496         /*Fetch number of nodes and dof for this finite element*/
     6497        /*Fetch number of nodes for this finite element*/
    64976498        int numnodes = this->NumberofNodes();
    6498         int numdof   = numnodes*NDOF1;
    64996499
    65006500        /*Initialize Element matrix and vectors*/
    65016501        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    65026502        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    6503         IssmDouble*    B      = xNew<IssmDouble>(2*numdof);
    6504         IssmDouble*    Bprime = xNew<IssmDouble>(2*numdof);
     6503        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
     6504        IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
    65056505        IssmDouble     D[2][2];
    65066506
    … …  
    65386538                D_scalar=gauss->weight*Jdettria;
    65396539
    6540                 TripleMultiply(basis,1,numdof,1,
     6540                TripleMultiply(basis,1,numnodes,1,
    65416541                                        &D_scalar,1,1,0,
    6542                                         basis,1,numdof,0,
     6542                                        basis,1,numnodes,0,
    65436543                                        &Ke->values[0],1);
    65446544
    … …  
    65546554                D[1][1]=D_scalar*dvydy;
    65556555                D[1][0]=0.;
    6556                 TripleMultiply(B,2,numdof,1,
     6556                TripleMultiply(B,2,numnodes,1,
    65576557                                        &D[0][0],2,2,0,
    6558                                         B,2,numdof,0,
     6558                                        B,2,numnodes,0,
    65596559                                        &Ke->values[0],1);
    65606560
    65616561                D[0][0]=D_scalar*vx;
    65626562                D[1][1]=D_scalar*vy;
    6563                 TripleMultiply(B,2,numdof,1,
     6563                TripleMultiply(B,2,numnodes,1,
    65646564                                        &D[0][0],2,2,0,
    6565                                         Bprime,2,numdof,0,
     6565                                        Bprime,2,numnodes,0,
    65666566                                        &Ke->values[0],1);
    65676567
    … …  
    65886588                        D[0][1]=D_scalar*D[0][1];
    65896589                        D[1][1]=D_scalar*D[1][1];
    6590                         TripleMultiply(Bprime,2,numdof,1,
     6590                        TripleMultiply(Bprime,2,numnodes,1,
    65916591                                                &D[0][0],2,2,0,
    6592                                                 Bprime,2,numdof,0,
     6592                                                Bprime,2,numnodes,0,
    65936593                                                &Ke->values[0],1);
    65946594                }
    … …  
    66076607
    66086608        /*Constants*/
    6609         const int    numdof=NDOF1*NUMVERTICES;
     6609        const int    numnodes=NUMVERTICES;
    66106610
    66116611        /*Intermediaries */
    … …  
    66536653                DL_scalar=gauss->weight*Jdettria;
    66546654
    6655                 TripleMultiply(&basis[0],1,numdof,1,
     6655                TripleMultiply(&basis[0],1,numnodes,1,
    66566656                                        &DL_scalar,1,1,0,
    6657                                         &basis[0],1,numdof,0,
     6657                                        &basis[0],1,numnodes,0,
    66586658                                        &Ke->values[0],1);
    66596659
    … …  
    66676667                DLprime[1][1]=DL_scalar*vy;
    66686668
    6669                 TripleMultiply( &B[0][0],2,numdof,1,
     6669                TripleMultiply( &B[0][0],2,numnodes,1,
    66706670                                        &DLprime[0][0],2,2,0,
    6671                                         &Bprime[0][0],2,numdof,0,
     6671                                        &Bprime[0][0],2,numnodes,0,
    66726672                                        &Ke->values[0],1);
    66736673        }
    … …  
    67016701        /*Fetch number of nodes and dof for this finite element*/
    67026702        int numnodes = this->NumberofNodes();
    6703         int numdof   = numnodes*NDOF1;
    67046703
    67056704        /*Initialize Element vector and other vectors*/
    … …  
    67336732                 basal_melting_correction_g=0.;
    67346733
    6735                 for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*basis[i];
     6734                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*basis[i];
    67366735        }
    67376736
    … …  
    67466745
    67476746        /*Constants*/
    6748         const int    numdof=NDOF1*NUMVERTICES;
     6747        const int    numnodes=NUMVERTICES;
    67496748
    67506749        /*Intermediaries */
    … …  
    67786777                thickness_input->GetInputValue(&thickness_g,gauss);
    67796778
    6780                 for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*basis[i];
     6779                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*basis[i];
    67816780        }
    67826781
    … …  
    69576956        /*Fetch number of nodes and dof for this finite element*/
    69586957        int numnodes = this->NumberofNodes();
    6959         int numdof   = numnodes*NDOF1;
    69606958
    69616959        /*Initialize Element matrix and vectors*/
    69626960        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
    69636961        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    6964         IssmDouble*    B      = xNew<IssmDouble>(2*numdof);
    6965         IssmDouble*    Bprime = xNew<IssmDouble>(2*numdof);
     6962        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
     6963        IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
    69666964        IssmDouble     D[2][2];
    69676965
    … …  
    70057003                D[1][1]=D_scalar*dvydy;
    70067004                D[1][1]=0.;
    7007                 TripleMultiply(B,2,numdof,1,
     7005                TripleMultiply(B,2,numnodes,1,
    70087006                                        &D[0][0],2,2,0,
    7009                                         B,2,numdof,0,
     7007                                        B,2,numnodes,0,
    70107008                                        &Ke->values[0],1);
    70117009
    70127010                D[0][0]=D_scalar*vx;
    70137011                D[1][1]=D_scalar*vy;
    7014                 TripleMultiply(B,2,numdof,1,
     7012                TripleMultiply(B,2,numnodes,1,
    70157013                                        &D[0][0],2,2,0,
    7016                                         Bprime,2,numdof,0,
     7014                                        Bprime,2,numnodes,0,
    70177015                                        &Ke->values[0],1);
    70187016
    … …  
    70397037                        D[0][1]=D_scalar*D[0][1];
    70407038                        D[1][1]=D_scalar*D[1][1];
    7041                         TripleMultiply(Bprime,2,numdof,1,
     7039                        TripleMultiply(Bprime,2,numnodes,1,
    70427040                                                &D[0][0],2,2,0,
    7043                                                 Bprime,2,numdof,0,
     7041                                                Bprime,2,numnodes,0,
    70447042                                                &Ke->values[0],1);
    70457043                }
    … …  
    70587056
    70597057        /*Constants*/
    7060         const int  numdof=NDOF1*NUMVERTICES;
     7058        const int  numnodes=NUMVERTICES;
    70617059
    70627060        /*Intermediaries*/
    … …  
    70977095                DL[1][1]=DL_scalar*vy;
    70987096
    7099                 TripleMultiply( &B[0][0],2,numdof,1,
     7097                TripleMultiply( &B[0][0],2,numnodes,1,
    71007098                                        &DL[0][0],2,2,0,
    7101                                         &Bprime[0][0],2,numdof,0,
     7099                                        &Bprime[0][0],2,numnodes,0,
    71027100                                        &Ke->values[0],1);
    71037101        }
    … …  
    71317129        /*Fetch number of nodes and dof for this finite element*/
    71327130        int numnodes = this->NumberofNodes();
    7133         int numdof   = numnodes*NDOF1;
    71347131
    71357132        /*Initialize Element vector and other vectors*/
    … …  
    71567153                GetNodalFunctions(basis,gauss);
    71577154
    7158                 for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
     7155                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
    71597156        }
    71607157
    … …  
    71697166
    71707167        /*Constants*/
    7171         const int    numdof=NDOF1*NUMVERTICES;
     7168        const int    numnodes=NUMVERTICES;
    71727169
    71737170        /*Intermediaries */
    … …  
    72007197                GetNodalFunctions(&basis[0],gauss);
    72017198
    7202                 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
     7199                for(i=0;i<numnodes;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
    72037200        }
    72047201
Note: See TracChangeset for help on using the changeset viewer.