Changeset 15730


Ignore:
Timestamp:
08/06/13 15:26:56 (12 years ago)
Author:
seroussi
Message:

FIX: fixing tiling HOFS, KMatrix first

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r15726 r15730  
    954954        /*First, figure out size of doflist and create it: */
    955955        int numberofdofs=0;
    956         for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
     956        for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum);
    957957
    958958        /*Allocate output*/
    … …  
    962962        int count=0;
    963963        for(int i=0;i<numnodes;i++){
    964                 nodes[i]->GetDofList(doflist+count,FSApproximationEnum,setenum);
    965                 count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
     964                nodes[i]->GetDofList(doflist+count,FSvelocityEnum,setenum);
     965                count+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum);
    966966        }
    967967
    … …  
    70907090        int       cs_list[numnodes];
    70917091        int       i,j;
     7092        int       init;
    70927093
    70937094        /*Prepare node list*/
    … …  
    70997100        }
    71007101
     7102        /*Some parameters needed*/
     7103        init=this->element_type;
     7104
    71017105        /*compute all stiffness matrices for this element*/
    71027106        ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
    … …  
    71057109        delete Ke1;
    71067110        delete Ke2;
    7107         Ke1=CreateKMatrixDiagnosticHO(); TransformInvStiffnessMatrixCoord(Ke1,this->nodes,NUMVERTICES,XYEnum);
    7108         Ke2=CreateKMatrixDiagnosticFS(); TransformInvStiffnessMatrixCoord(Ke2,this->nodes,NUMVERTICES,XYZEnum);
     7111       
     7112        this->element_type=P1Enum;
     7113        Ke1=CreateKMatrixDiagnosticHO();
     7114        this->element_type=init;
     7115        Ke2=CreateKMatrixDiagnosticFS();
    71097116
    71107117        for(i=0;i<numdofs;i++) for(j=0;j<NUMVERTICES;j++){
    … …  
    71317138        int approximation;
    71327139        inputs->GetInputValue(&approximation,ApproximationEnum);
    7133 
    71347140        switch(approximation){
    71357141                case SSAApproximationEnum:
    … …  
    76647670
    76657671        /*compute all stiffness matrices for this element*/
     7672        int init = this->element_type;
     7673        this->element_type=P1Enum;
    76667674        ElementMatrix* Ke1=CreateKMatrixDiagnosticHO();
     7675        this->element_type=init;
    76677676        ElementMatrix* Ke2=CreateKMatrixDiagnosticFS();
    76687677        ElementMatrix* Ke3=CreateKMatrixCouplingHOFS();
    … …  
    95309539        for(int i=0;i<vnumnodes;i++){
    95319540                gauss->GaussNode(this->VelocityInterpolation(),i);
    9532 
    95339541                vx_input->GetInputValue(&vx,gauss);
    95349542                vy_input->GetInputValue(&vy,gauss);
    … …  
    96709678void  Penta::InputUpdateFromSolutionDiagnosticSSA(IssmDouble* solution){
    96719679
    9672         const int    numdof=NDOF2*NUMVERTICES;
    9673 
    9674         int     i;
     9680        int         numnodes = this->NumberofNodes();
     9681        int         numdof=NDOF2*numnodes;
     9682
     9683        int         i;
    96759684        IssmDouble  rho_ice,g;
    96769685        IssmDouble  values[numdof];
    … …  
    98619870        /*Get dof listof this element (SSA dofs) and of the penta at base (SSA dofs): */
    98629871        penta->GetDofList(&doflistm,SSAApproximationEnum,GsetEnum);
    9863         GetDofList(&doflists,FSApproximationEnum,GsetEnum);
     9872        GetDofList(&doflists,FSvelocityEnum,GsetEnum);
    98649873        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    98659874
    … …  
    1009010099
    1009110100        const int    numdofp=NDOF2*NUMVERTICES;
    10092         const int    numdofs=NDOF4*NUMVERTICES;
    10093 
    10094         int    i;
     10101        const int    numdofs=NDOF3*NUMVERTICES;
     10102        const int    numdofpressure=NDOF1*NUMVERTICES;
     10103
     10104        int        i;
    1009510105        IssmDouble HO_values[numdofp];
    1009610106        IssmDouble FS_values[numdofs];
     10107        IssmDouble Pressure_values[numdofpressure];
    1009710108        IssmDouble vx[NUMVERTICES];
    1009810109        IssmDouble vy[NUMVERTICES];
    … …  
    1010410115        IssmDouble xyz_list[NUMVERTICES][3];
    1010510116        IssmDouble FSreconditioning;
    10106         int*   doflistp      = NULL;
    10107         int*   doflists      = NULL;
    10108         Penta  *penta        = NULL;
     10117        int*       doflistp        = NULL;
     10118        int*       doflists        = NULL;
     10119        int*       doflistpressure = NULL;
     10120        Penta      *penta          = NULL;
    1010910121
    1011010122        /*OK, we have to add results of this element for HO
    … …  
    1011410126        /*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
    1011510127        GetDofList(&doflistp,HOApproximationEnum,GsetEnum);
    10116         GetDofList(&doflists,FSApproximationEnum,GsetEnum);
     10128        GetDofList(&doflists,FSvelocityEnum,GsetEnum);
     10129        GetDofListPressure(&doflistpressure,GsetEnum);
    1011710130        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    1011810131
    … …  
    1012310136        for(i=0;i<numdofp;i++) HO_values[i]=solution[doflistp[i]];
    1012410137        for(i=0;i<numdofs;i++) FS_values[i]=solution[doflists[i]];
     10138        for(i=0;i<numdofpressure;i++) Pressure_values[i]=solution[doflistpressure[i]];
    1012510139
    1012610140        /*Transform solution in Cartesian Space*/
    … …  
    1013010144        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    1013110145        for(i=0;i<NUMVERTICES;i++){
    10132                 vx[i]=FS_values[i*NDOF4+0]+HO_values[i*NDOF2+0];
    10133                 vy[i]=FS_values[i*NDOF4+1]+HO_values[i*NDOF2+1];
    10134                 vzFS[i]=FS_values[i*NDOF4+2];
    10135                 pressure[i]=FS_values[i*NDOF4+3]*FSreconditioning;
     10146                vx[i]=FS_values[i*NDOF3+0]+HO_values[i*NDOF2+0];
     10147                vy[i]=FS_values[i*NDOF3+1]+HO_values[i*NDOF2+1];
     10148                vzFS[i]=FS_values[i*NDOF3+2];
     10149                pressure[i]=FS_values[i*NDOF1]*FSreconditioning;
    1013610150
    1013710151                /*Check solution*/
    … …  
    1017810192        xDelete<int>(doflistp);
    1017910193        xDelete<int>(doflists);
     10194        xDelete<int>(doflistpressure);
    1018010195}
    1018110196/*}}}*/
    … …  
    1018310198void  Penta::InputUpdateFromSolutionDiagnosticSIA(IssmDouble* solution){
    1018410199
    10185         const int    numdof=NDOF2*NUMVERTICES;
     10200        int         numnodes = this->NumberofNodes();
     10201        int         numdof=NDOF2*numnodes;
    1018610202
    1018710203        int     i;
    … …  
    1024610262void  Penta::InputUpdateFromSolutionDiagnosticVert(IssmDouble* solution){
    1024710263
    10248         const int numdof=NDOF1*NUMVERTICES;
    10249 
    10250         int      i;
    10251         int      approximation;
     10264        int          numnodes = this->NumberofNodes();
     10265        int          numdof=NDOF1*numnodes;
     10266
     10267        int          i;
     10268        int          approximation;
    1025210269        IssmDouble   rho_ice,g;
    1025310270        IssmDouble   values[numdof];
    … …  
    1026210279        IssmDouble   surface[NUMVERTICES];
    1026310280        IssmDouble   xyz_list[NUMVERTICES][3];
    10264         int*     doflist      = NULL;
     10281        int*         doflist      = NULL;
    1026510282
    1026610283        /*Get the approximation and do nothing if the element in FS or None*/
Note: See TracChangeset for help on using the changeset viewer.