Changeset 16433


Ignore:
Timestamp:
10/16/13 13:34:51 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: fixed coupling FS

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

Legend:

Unmodified
Added
Removed

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

    r16422 r16433  
    1016010160void  Penta::InputUpdateFromSolutionStressbalanceSSAFS(IssmDouble* solution){
    1016110161
    10162         const int    numdofm=NDOF2*NUMVERTICES;
    10163         const int    numdofs=NDOF3*NUMVERTICES;
    10164         const int    numdof2d=NDOF2*NUMVERTICES2D;
    10165         const int    numdofpressure=NDOF1*NUMVERTICES;
     10162        const int    numdofSSA = NDOF2*NUMVERTICES;
     10163        const int    numdof2d  = NDOF2*NUMVERTICES2D;
     10164        const int    numdofFSv = NDOF3*NUMVERTICES;
     10165        const int    numdofFSp = NDOF1*NUMVERTICES;
    1016610166
    1016710167        int     i;
    1016810168        IssmDouble  FSreconditioning;
    10169         IssmDouble  SSA_values[numdofm];
    10170         IssmDouble  FS_values[numdofs];
    10171         IssmDouble  Pressure_values[numdofs];
     10169        IssmDouble  SSA_values[numdofSSA];
     10170        IssmDouble  FS_values[numdofFSv+numdofFSp];
    1017210171        IssmDouble  vx[NUMVERTICES];
    1017310172        IssmDouble  vy[NUMVERTICES];
    … …  
    1017810177        IssmDouble  pressure[NUMVERTICES];
    1017910178        IssmDouble  xyz_list[NUMVERTICES][3];
    10180         int*    doflistm        = NULL;
    10181         int*    doflists        = NULL;
    10182         int*    doflistpressure = NULL;
    10183         Penta   *penta          = NULL;
     10179        int   *doflistSSA = NULL;
     10180        int   *doflistFSv = NULL;
     10181        int   *doflistFSp = NULL;
     10182        Penta *penta      = NULL;
     10183
     10184        /*Prepare coordinate system list*/
     10185        int* cs_list = xNew<int>(2*NUMVERTICES);
     10186        for(i=0;i<NUMVERTICES;i++) cs_list[i]             = XYZEnum;
     10187        for(i=0;i<NUMVERTICES;i++) cs_list[NUMVERTICES+i] = PressureEnum;
    1018410188
    1018510189        /*OK, we have to add results of this element for SSA
    … …  
    1018810192
    1018910193        /*Get dof listof this element (SSA dofs) and of the penta at base (SSA dofs): */
    10190         penta->GetDofList(&doflistm,SSAApproximationEnum,GsetEnum);
    10191         GetDofList(&doflists,FSvelocityEnum,GsetEnum);
    10192         GetDofListPressure(&doflistpressure,GsetEnum);
     10194        penta->GetDofList(&doflistSSA,SSAApproximationEnum,GsetEnum);
     10195        GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
     10196        GetDofListPressure(&doflistFSp,GsetEnum);
    1019310197        this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    1019410198
    … …  
    1019810202        /*Use the dof list to index into the solution vector: */
    1019910203        for(i=0;i<numdof2d;i++){
    10200                 SSA_values[i]=solution[doflistm[i]];
    10201                 SSA_values[i+numdof2d]=solution[doflistm[i]];
    10202         }
    10203         for(i=0;i<numdofs;i++)FS_values[i]=solution[doflists[i]];
    10204         for(i=0;i<numdofpressure;i++) Pressure_values[i]=solution[doflistpressure[i]];
     10204                SSA_values[i]=solution[doflistSSA[i]];
     10205                SSA_values[i+numdof2d]=solution[doflistSSA[i]];
     10206        }
     10207        for(i=0;i<numdofFSv;i++) FS_values[i]=solution[doflistFSv[i]];
     10208        for(i=0;i<numdofFSp;i++) FS_values[numdofFSv+i]=solution[doflistFSp[i]];
    1020510209
    1020610210        /*Transform solution in Cartesian Space*/
    1020710211        TransformSolutionCoord(&SSA_values[0],this->nodes,NUMVERTICES,XYEnum);
    10208         TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZEnum);
     10212        TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
    1020910213
    1021010214        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    1021110215        for(i=0;i<NUMVERTICES;i++){
    10212                 vx[i]=FS_values[i*NDOF3+0]+SSA_values[i*NDOF2+0];
    10213                 vy[i]=FS_values[i*NDOF3+1]+SSA_values[i*NDOF2+1];
    10214                 vzFS[i]=FS_values[i*NDOF3+2];
    10215                 pressure[i]=Pressure_values[i*NDOF1]*FSreconditioning;
     10216                vx[i]       = FS_values[i*NDOF3+0]+SSA_values[i*NDOF2+0];
     10217                vy[i]       = FS_values[i*NDOF3+1]+SSA_values[i*NDOF2+1];
     10218                vzFS[i]     = FS_values[i*NDOF3+2];
     10219                pressure[i] = FS_values[NUMVERTICES*NDOF3+i]*FSreconditioning;
    1021610220
    1021710221                /*Check solution*/
    1021810222                if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
    1021910223                if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
    10220                 if(xIsNan<IssmDouble>(vzFS[i])) _error_("NaN found in solution vector");
     10224                if(xIsNan<IssmDouble>(vzFS[i]))     _error_("NaN found in solution vector");
    1022110225                if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
    1022210226        }
    … …  
    1023610240        /*Now Compute vel*/
    1023710241        for(i=0;i<NUMVERTICES;i++) {
    10238                 vz[i]=vzSSA[i]+vzFS[i];
    10239                 vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
     10242                vz[i]  = vzSSA[i]+vzFS[i];
     10243                vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    1024010244        }
    1024110245
    … …  
    1025610260
    1025710261        /*Free ressources:*/
    10258         xDelete<int>(doflistm);
    10259         xDelete<int>(doflists);
    10260         xDelete<int>(doflistpressure);
     10262        xDelete<int>(doflistSSA);
     10263        xDelete<int>(doflistFSv);
     10264        xDelete<int>(doflistFSp);
     10265        xDelete<int>(cs_list);
    1026110266}
    1026210267/*}}}*/
    … …  
    1041810423void  Penta::InputUpdateFromSolutionStressbalanceHOFS(IssmDouble* solution){
    1041910424
    10420         const int    numdofp=NDOF2*NUMVERTICES;
    10421         const int    numdofs=NDOF3*NUMVERTICES;
    10422         const int    numdofpressure=NDOF1*NUMVERTICES;
     10425        const int    numdofHO  = NDOF2*NUMVERTICES;
     10426        const int    numdofFSv = NDOF3*NUMVERTICES;
     10427        const int    numdofFSp = NDOF1*NUMVERTICES;
    1042310428
    1042410429        int        i;
    10425         IssmDouble HO_values[numdofp];
    10426         IssmDouble FS_values[numdofs];
    10427         IssmDouble Pressure_values[numdofpressure];
     10430        IssmDouble HO_values[numdofHO];
     10431        IssmDouble FS_values[numdofFSv+numdofFSp];
    1042810432        IssmDouble vx[NUMVERTICES];
    1042910433        IssmDouble vy[NUMVERTICES];
    … …  
    1043510439        IssmDouble xyz_list[NUMVERTICES][3];
    1043610440        IssmDouble FSreconditioning;
    10437         int*       doflistp        = NULL;
    10438         int*       doflists        = NULL;
    10439         int*       doflistpressure = NULL;
     10441        int*       doflistHO        = NULL;
     10442        int*       doflistFSv        = NULL;
     10443        int*       doflistFSp = NULL;
     10444
     10445        /*Prepare coordinate system list*/
     10446        int* cs_list = xNew<int>(2*NUMVERTICES);
     10447        for(i=0;i<NUMVERTICES;i++) cs_list[i]             = XYZEnum;
     10448        for(i=0;i<NUMVERTICES;i++) cs_list[NUMVERTICES+i] = PressureEnum;
    1044010449
    1044110450        /*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
    10442         GetDofList(&doflistp,HOApproximationEnum,GsetEnum);
    10443         GetDofList(&doflists,FSvelocityEnum,GsetEnum);
    10444         GetDofListPressure(&doflistpressure,GsetEnum);
     10451        GetDofList(&doflistHO,HOApproximationEnum,GsetEnum);
     10452        GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
     10453        GetDofListPressure(&doflistFSp,GsetEnum);
    1044510454        this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    1044610455
    … …  
    1044910458
    1045010459        /*Use the dof list to index into the solution vector: */
    10451         for(i=0;i<numdofp;i++) HO_values[i]=solution[doflistp[i]];
    10452         for(i=0;i<numdofs;i++) FS_values[i]=solution[doflists[i]];
    10453         for(i=0;i<numdofpressure;i++) Pressure_values[i]=solution[doflistpressure[i]];
     10460        for(i=0;i<numdofHO;i++)  HO_values[i]=solution[doflistHO[i]];
     10461        for(i=0;i<numdofFSv;i++) FS_values[i]=solution[doflistFSv[i]];
     10462        for(i=0;i<numdofFSp;i++) FS_values[numdofFSv+i]=solution[doflistFSp[i]];
    1045410463
    1045510464        /*Transform solution in Cartesian Space*/
    1045610465        TransformSolutionCoord(&HO_values[0],this->nodes,NUMVERTICES,XYEnum);
    10457         TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZEnum);
     10466        TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
    1045810467
    1045910468        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    1046010469        for(i=0;i<NUMVERTICES;i++){
    10461                 vx[i]=FS_values[i*NDOF3+0]+HO_values[i*NDOF2+0];
    10462                 vy[i]=FS_values[i*NDOF3+1]+HO_values[i*NDOF2+1];
    10463                 vzFS[i]=FS_values[i*NDOF3+2];
    10464                 pressure[i]=Pressure_values[i*NDOF1]*FSreconditioning;
     10470                vx[i]       = FS_values[i*NDOF3+0]+HO_values[i*NDOF2+0];
     10471                vy[i]       = FS_values[i*NDOF3+1]+HO_values[i*NDOF2+1];
     10472                vzFS[i]     = FS_values[i*NDOF3+2];
     10473                pressure[i] = FS_values[NUMVERTICES*NDOF3+i]*FSreconditioning;
    1046510474
    1046610475                /*Check solution*/
    1046710476                if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
    1046810477                if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
    10469                 if(xIsNan<IssmDouble>(vzFS[i])) _error_("NaN found in solution vector");
     10478                if(xIsNan<IssmDouble>(vzFS[i]))     _error_("NaN found in solution vector");
    1047010479                if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
    1047110480        }
    … …  
    1048510494        /*Now Compute vel*/
    1048610495        for(i=0;i<NUMVERTICES;i++) {
    10487                 vz[i]=vzHO[i]+vzFS[i];
    10488                 vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
     10496                vz[i]  = vzHO[i]+vzFS[i];
     10497                vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    1048910498        }
    1049010499
    … …  
    1050510514
    1050610515        /*Free ressources:*/
    10507         xDelete<int>(doflistp);
    10508         xDelete<int>(doflists);
    10509         xDelete<int>(doflistpressure);
     10516        xDelete<int>(doflistHO);
     10517        xDelete<int>(doflistFSv);
     10518        xDelete<int>(doflistFSp);
     10519        xDelete<int>(cs_list);
    1051010520}
    1051110521/*}}}*/
    … …  
    1070210712        /*Prepare coordinate system list*/
    1070310713        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    10704         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
     10714        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
    1070510715        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    1070610716

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

    r16429 r16433  
    24552455        GaussTria* gauss = new GaussTria();
    24562456        for(int i=0;i<numindices;i++){//FIXME
    2457 
    24582457
    24592458                gauss->GaussNode(this->VelocityInterpolation(),indices[i]);
Note: See TracChangeset for help on using the changeset viewer.