Ignore:
Timestamp:
05/05/09 14:11:46 (17 years ago)
Author:
Eric.Larour
Message:

New parameter inputs framework, that can adapt to unknow number of dofs problems

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r241 r246  
    266266#define __FUNCT__ "Tria::CreateKMatrix"
    267267
    268 void  Tria::CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type){
     268void  Tria::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){
    269269
    270270        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    … …  
    284284#define __FUNCT__ "Tria::CreateKMatrixDiagnosticHoriz"
    285285
    286 void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,ParameterInputs* inputs,int analysis_type){
     286void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type){
    287287
    288288
    … …  
    328328       
    329329        /*input parameters for structural analysis (diagnostic): */
    330         double* velocity_param=NULL;
    331         double* oldvelocity_param=NULL;
    332         double  vxvy_list[numgrids][2];
    333         double  oldvxvy_list[numgrids][2];
     330        double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
     331        double  oldvxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
    334332        double  thickness;
     333        int     dofs[2]={0,1};
     334
     335        ParameterInputs* inputs=NULL;
     336
     337        /*recover pointers: */
     338        inputs=(ParameterInputs*)vinputs;
    335339
    336340        /*recover extra inputs from users, at current convergence iteration: */
    337         velocity_param=ParameterInputsRecover(inputs,"velocity");
    338         oldvelocity_param=ParameterInputsRecover(inputs,"old_velocity");
     341        inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
     342        inputs->Recover("old_velocity",&oldvxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
    339343
    340344        /* Get node coordinates and dof list: */
    … …  
    344348        /* Set Ke_gg to 0: */
    345349        for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
    346 
    347         /*Initialize vxvy_list: */
    348         for(i=0;i<numgrids;i++){
    349                 if(velocity_param){
    350                         vxvy_list[i][0]=velocity_param[doflist[i*numberofdofspernode+0]];
    351                         vxvy_list[i][1]=velocity_param[doflist[i*numberofdofspernode+1]];
    352                 }
    353                 else{
    354                         vxvy_list[i][0]=0;
    355                         vxvy_list[i][1]=0;
    356                 }
    357 
    358                 if(oldvelocity_param){
    359                         oldvxvy_list[i][0]=oldvelocity_param[doflist[i*numberofdofspernode+0]];
    360                         oldvxvy_list[i][1]=oldvelocity_param[doflist[i*numberofdofspernode+1]];
    361                 }
    362                 else{
    363                         oldvxvy_list[i][0]=0;
    364                         oldvxvy_list[i][1]=0;
    365                 }
    366         }
    367350
    368351        #ifdef _DEBUGELEMENTS_
    … …  
    516499#undef __FUNCT__
    517500#define __FUNCT__ "Tria::CreateKMatrixDiagnosticHorizFriction"
    518 void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,ParameterInputs* inputs,int analysis_type){
     501void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type){
    519502
    520503
    … …  
    554537
    555538        /*input parameters for structural analysis (diagnostic): */
    556         double* velocity_param=NULL;
    557         double  vxvy_list[numgrids][2];
     539        double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
     540        int     dofs[2]={0,1};
    558541
    559542        /*friction: */
    … …  
    564547        double MOUNTAINKEXPONENT=10;
    565548
    566         /*recover extra inputs from users, at current convergence iteration: */
    567         velocity_param=ParameterInputsRecover(inputs,"velocity");
    568 
     549        ParameterInputs* inputs=NULL;
     550
     551        /*recover pointers: */
     552        inputs=(ParameterInputs*)vinputs;
     553       
    569554        /* Get node coordinates and dof list: */
    570555        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    … …  
    581566        if (friction_type!=2)throw ErrorException(__FUNCT__," non-viscous friction not supported yet!");
    582567
    583         /*Initialize vxvy_list: */
    584         for(i=0;i<numgrids;i++){
    585                 if(velocity_param){
    586                         vxvy_list[i][0]=velocity_param[doflist[i*numberofdofspernode+0]];
    587                         vxvy_list[i][1]=velocity_param[doflist[i*numberofdofspernode+1]];
    588                 }
    589                 else{
    590                         vxvy_list[i][0]=0;
    591                         vxvy_list[i][1]=0;
    592                 }
    593         }
     568        /*recover extra inputs from users, at current convergence iteration: */
     569        inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
    594570
    595571        /*Build alpha2_list used by drag stiffness matrix*/
    … …  
    691667#undef __FUNCT__
    692668#define __FUNCT__ "Tria::CreatePVector"
    693 void  Tria::CreatePVector(Vec pg,ParameterInputs* inputs,int analysis_type){
     669void  Tria::CreatePVector(Vec pg,void* inputs,int analysis_type){
    694670       
    695671        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
    … …  
    708684#define __FUNCT__ "Tria::CreatePVectorDiagnosticHoriz"
    709685
    710 void Tria::CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs, int analysis_type){
     686void Tria::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs, int analysis_type){
    711687
    712688        int             i,j;
    … …  
    746722        /*input parameters for structural analysis (diagnostic): */
    747723        double  thickness;
     724
     725        ParameterInputs* inputs=NULL;
     726
     727        /*recover pointers: */
     728        inputs=(ParameterInputs*)vinputs;
    748729
    749730        /* Get node coordinates and dof list: */
    … …  
    871852#undef __FUNCT__
    872853#define __FUNCT__ "Tria::UpdateFromInputs"
    873 void  Tria::UpdateFromInputs(ParameterInputs* inputs){
    874 
    875         int i;
    876         int doflist[MAXDOFSPERNODE*3];
    877         int numberofdofspernode;
    878 
    879         /*element: */
    880         double* h_param=NULL;
    881         double* s_param=NULL;
    882         double* b_param=NULL;
    883         double* k_param=NULL;
    884         double* melting_param=NULL;
    885         double* accumulation_param=NULL;
    886        
    887         /*material: */
    888         double* temperature_average_param=NULL;
     854void  Tria::UpdateFromInputs(void* vinputs){
     855
     856        int     dofs[1]={0};
     857        double  temperature_list[3];
    889858        double  temperature_average;
    890         double* B_param=NULL;
     859        double  B_list[3];
    891860        double  B_average;
    892861
    893 
    894         /*Get dof list: */
    895         GetDofList(&doflist[0],&numberofdofspernode);
     862        ParameterInputs* inputs=NULL;
     863
     864        /*recover pointers: */
     865        inputs=(ParameterInputs*)vinputs;
    896866
    897867        /*Update internal data if inputs holds new values: */
    898         h_param=ParameterInputsRecover(inputs,"thickness");
    899         s_param=ParameterInputsRecover(inputs,"surface");
    900         b_param=ParameterInputsRecover(inputs,"bed");
    901         k_param=ParameterInputsRecover(inputs,"drag");
    902         melting_param=ParameterInputsRecover(inputs,"melting");
    903         accumulation_param=ParameterInputsRecover(inputs,"accumulation");
    904 
    905         for(i=0;i<3;i++){
    906                 if(h_param)h[i]=h_param[doflist[i*numberofdofspernode+0]];
    907                 if(s_param)s[i]=s_param[doflist[i*numberofdofspernode+0]];
    908                 if(b_param)b[i]=b_param[doflist[i*numberofdofspernode+0]];
    909                 if(k_param)k[i]=k_param[doflist[i*numberofdofspernode+0]];
    910                 if(melting_param)melting[i]=melting_param[doflist[i*numberofdofspernode+0]];
    911                 if(accumulation_param)accumulation[i]=accumulation_param[doflist[i*numberofdofspernode+0]];
    912         }
    913 
     868        inputs->Recover("thickness",&h[0],1,dofs,3,(void**)nodes);
     869        inputs->Recover("surface",&s[0],1,dofs,3,(void**)nodes);
     870        inputs->Recover("bed",&b[0],1,dofs,3,(void**)nodes);
     871        inputs->Recover("drag",&k[0],1,dofs,3,(void**)nodes);
     872        inputs->Recover("melting",&melting[0],1,dofs,3,(void**)nodes);
     873        inputs->Recover("accumulation_param",&accumulation[0],1,dofs,3,(void**)nodes);
     874       
    914875        //Update material if necessary
    915         temperature_average_param=ParameterInputsRecover(inputs,"temperature_average");
    916         B_param=ParameterInputsRecover(inputs,"B");
    917         if(temperature_average_param){
    918                 temperature_average-0;
    919                 for(i=0;i<3;i++) temperature_average+=temperature_average_param[doflist[i*numberofdofspernode+0]];
    920                 temperature_average=temperature_average/3.0;
     876        if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,3,(void**)nodes)){
     877                temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2])/3.0;
    921878                B_average=Paterson(temperature_average);
    922879                matice->SetB(B_average);
    923880        }
    924 
    925         if(B_param){
    926                 B_average=0;
    927                 for(i=0;i<3;i++) B_average+=B_param[doflist[i*numberofdofspernode+0]];
    928                 B_average=B_average/3.0;
     881       
     882        if(inputs->Recover("B",&B_list[0],1,dofs,3,(void**)nodes)){
     883                B_average=(B_list[0]+B_list[1]+B_list[2])/3.0;
    929884                matice->SetB(B_average);
    930885        }
    … …  
    13111266}
    13121267
    1313 Matpar* Tria::GetMatPar(){
     1268void* Tria::GetMatPar(){
    13141269        return matpar;
    13151270}
    … …  
    13201275
    13211276               
    1322 void  Tria::GetNodes(Node** pnodes){
     1277void  Tria::GetNodes(void** vpnodes){
    13231278        int i;
     1279        Node** pnodes=(Node**)vpnodes;
     1280
    13241281        for(i=0;i<3;i++){
    13251282                pnodes[i]=nodes[i];
    … …  
    13461303
    13471304Object* Tria::copy() {
     1305       
    13481306        return new Tria(*this);
     1307
    13491308}
    13501309
    … …  
    13521311#undef __FUNCT__
    13531312#define __FUNCT__ "Tria::Du"
    1354 void Tria::Du(Vec du_g,double* velocity,double* obs_velocity,ParameterInputs* inputs,int analysis_type){
     1313void Tria::Du(Vec du_g,double* velocity,double* obs_velocity,void* vinputs,int analysis_type){
    13551314
    13561315        int i;
    … …  
    14041363        double scaley=0;
    14051364        double scale=0;
    1406         double* fit_param=NULL;
    14071365        double  fit=-1;
     1366
     1367        ParameterInputs* inputs=NULL;
     1368
     1369        /*recover pointers: */
     1370        inputs=(ParameterInputs*)vinputs;
    14081371
    14091372        /* Get node coordinates and dof list: */
    … …  
    14151378
    14161379        /* Recover input data: */
    1417         fit_param=ParameterInputsRecover(inputs,"fit");
    1418         if(fit_param){
    1419                 fit=*fit_param;
    1420         }
    1421         else throw ErrorException(__FUNCT__," missing fit input parameter");
     1380        if(!inputs->Recover("fit",&fit)) throw ErrorException(__FUNCT__," missing fit input parameter");
    14221381
    14231382        /*Initialize velocities: */
    … …  
    15701529#undef __FUNCT__
    15711530#define __FUNCT__ "Tria::Gradj"
    1572 void  Tria::Gradj(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type,char* control_type){
     1531void  Tria::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,char* control_type){
    15731532
    15741533        if (strcmp(control_type,"drag")==0){
    … …  
    15831542#undef __FUNCT__
    15841543#define __FUNCT__ "Tria::GradjDrag"
    1585 void  Tria::GradjDrag(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type){
     1544void  Tria::GradjDrag(Vec grad_g,double* velocity,double* adjoint,void* vinputs,int analysis_type){
    15861545
    15871546
    … …  
    16031562        double adjy_list[numgrids];
    16041563
    1605         double* basal_drag_param=NULL;
    1606         double K_list[numgrids];
    16071564        double drag;
    1608 
    1609         double* bed_param=NULL;
    1610         double bed_list[numgrids];
    1611         double* thickness_param=NULL;
    1612         double thickness_list[numgrids];
     1565        int    dofs[1]={0};
    16131566
    16141567        /* gaussian points: */
    … …  
    16501603        double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    16511604
     1605        ParameterInputs* inputs=NULL;
     1606
     1607        /*recover pointers: */
     1608        inputs=(ParameterInputs*)vinputs;
     1609
    16521610        /* Get node coordinates and dof list: */
    16531611        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    … …  
    16581616
    16591617        /* recover input parameters: */
    1660         basal_drag_param=ParameterInputsRecover(inputs,"drag");
    1661         bed_param=ParameterInputsRecover(inputs,"bed");
    1662         thickness_param=ParameterInputsRecover(inputs,"thickness");
     1618        inputs->Recover("drag",&k[0],1,dofs,numgrids,(void**)nodes);
     1619        inputs->Recover("bed",&b[0],1,dofs,numgrids,(void**)nodes);
     1620        inputs->Recover("thickness",&h[0],1,dofs,numgrids,(void**)nodes);
    16631621
    16641622        /*Initialize parameter lists: */
    … …  
    16701628                adjx_list[i]=adjoint[doflist[i*numberofdofspernode+0]];
    16711629                adjy_list[i]=adjoint[doflist[i*numberofdofspernode+1]];
    1672                 if(basal_drag_param){
    1673                         K_list[i]=basal_drag_param[doflist[i*numberofdofspernode+0]];
    1674                 }
    1675                 else{
    1676                         K_list[i]=k[i];
    1677                 }
    1678                 if(bed_param){
    1679                         bed_list[i]=bed_param[doflist[i*numberofdofspernode+0]];
    1680                 }
    1681                 else{
    1682                         bed_list[i]=b[i];
    1683                 }
    1684                 if(thickness_param){
    1685                         thickness_list[i]=thickness_param[doflist[i*numberofdofspernode+0]];
    1686                 }
    1687                 else{
    1688                         thickness_list[i]=h[i];
    1689                 }
    16901630        }
    16911631       
    … …  
    17231663                        friction->rho_ice=matpar->GetRhoIce();
    17241664                        friction->rho_water=matpar->GetRhoWater();
    1725                         friction->K=&K_list[0];
    1726                         friction->bed=&bed_list[0];
    1727                         friction->thickness=&thickness_list[0];
     1665                        friction->K=&k[0];
     1666                        friction->bed=&b[0];
     1667                        friction->thickness=&h[0];
    17281668                        friction->velocities=&vxvy_list[0][0];
    17291669                        friction->p=p;
    … …  
    17461686                /*Recover alpha_complement and k: */
    17471687                GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3);
    1748                 GetParameterValue(&drag, &K_list[0],gauss_l1l2l3);
     1688                GetParameterValue(&drag, &k[0],gauss_l1l2l3);
    17491689                #ifdef _DEBUG_
    17501690                        printf("Drag complement: %20.20lf Drag: %20.20lf\n",alpha_complement,drag);
    … …  
    17961736#undef __FUNCT__
    17971737#define __FUNCT__ "Tria::GradjB"
    1798 void  Tria::GradjB(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type){
     1738void  Tria::GradjB(Vec grad_g,double* velocity,double* adjoint,void* vinputs,int analysis_type){
    17991739
    18001740        int i;
    … …  
    18021742        /* node data: */
    18031743        const int    numgrids=3;
    1804         const int    numdofs=2*numgrids;
    1805         const int    NDOF2=2;
     1744        const int    NDOF1=1;
     1745        const int    NDOF2=1;
     1746        const int    numdofs=NDOF1*numgrids;
    18061747        double       xyz_list[numgrids][3];
    18071748        int          doflist[numdofs];
    … …  
    18141755        double adjx_list[numgrids];
    18151756        double adjy_list[numgrids];
    1816 
    1817         double* thickness_param=NULL;
    1818         double thickness_list[numgrids];
    18191757
    18201758        /* gaussian points: */
    … …  
    18501788        double  lambda,mu;
    18511789        double  thickness;
    1852        
     1790        int     dofs[1]={0};
     1791       
     1792        ParameterInputs* inputs=NULL;
     1793
     1794        /*recover pointers: */
     1795        inputs=(ParameterInputs*)vinputs;
    18531796
    18541797        /* Get node coordinates and dof list: */
    … …  
    18601803
    18611804        /* recover input parameters: */
    1862         thickness_param=ParameterInputsRecover(inputs,"thickness");
     1805        inputs->Recover("thickness",&h[0],1,dofs,numgrids,(void**)nodes);
    18631806
    18641807        /*Initialize parameter lists: */
    … …  
    18701813                adjx_list[i]=adjoint[doflist[i*numberofdofspernode+0]];
    18711814                adjy_list[i]=adjoint[doflist[i*numberofdofspernode+1]];
    1872                 if(thickness_param){
    1873                         thickness_list[i]=thickness_param[doflist[i*numberofdofspernode+0]];
    1874                 }
    1875                 else{
    1876                         thickness_list[i]=h[i];
    1877                 }
    18781815        }
    18791816       
    … …  
    19401877#undef __FUNCT__
    19411878#define __FUNCT__ "Tria::Misfit"
    1942 double Tria::Misfit(double* velocity,double* obs_velocity,ParameterInputs* inputs,int analysis_type){
     1879double Tria::Misfit(double* velocity,double* obs_velocity,void* vinputs,int analysis_type){
    19431880        int i;
    19441881       
    … …  
    19791916        double scalex=1;
    19801917        double scaley=1;
    1981         double* fit_param=NULL;
    19821918        double  fit=-1;
     1919
     1920        ParameterInputs* inputs=NULL;
     1921
     1922        /*recover pointers: */
     1923        inputs=(ParameterInputs*)vinputs;
    19831924
    19841925        /* Get node coordinates and dof list: */
    … …  
    19871928       
    19881929        /* Recover input data: */
    1989         fit_param=ParameterInputsRecover(inputs,"fit");
    1990         if(fit_param){
    1991                 fit=*fit_param;
    1992         }
    1993         else ErrorException(__FUNCT__," missing fit input parameter");
     1930        if(!inputs->Recover("fit",&fit)) throw ErrorException(__FUNCT__," missing fit input parameter");
    19941931
    19951932        /*Initialize velocities: */
    … …  
    21022039
    21032040
     2041
     2042#undef __FUNCT__
     2043#define __FUNCT__ "Tria::CreateKMatrixDiagnosticBaseVert"
     2044void  Tria::CreateKMatrixDiagnosticBaseVert(Mat Kgg,void* vinputs,int analysis_type){
     2045
     2046        int i,j;
     2047
     2048        /* node data: */
     2049        const int    numgrids=3;
     2050        const int    NDOF1=1;
     2051        const int    numdofs=NDOF1*numgrids;
     2052        double       xyz_list[numgrids][3];
     2053        int          doflist[numdofs];
     2054        int          numberofdofspernode;
     2055       
     2056        /* gaussian points: */
     2057        int     num_gauss,ig;
     2058        double* first_gauss_area_coord  =  NULL;
     2059        double* second_gauss_area_coord =  NULL;
     2060        double* third_gauss_area_coord  =  NULL;
     2061        double* gauss_weights           =  NULL;
     2062        double  gauss_weight;
     2063        double  gauss_l1l2l3[3];
     2064
     2065
     2066        /* matrices: */
     2067        double L[1][3];
     2068        double DL_scalar;
     2069
     2070        double Ke_gg[3][3]; //stiffness matrix
     2071        double Ke_gg_gaussian[3][3]; //stiffness matrix evaluated at the gaussian point.
     2072        double Jdet;
     2073
     2074        ParameterInputs* inputs=NULL;
     2075
     2076        /*recover pointers: */
     2077        inputs=(ParameterInputs*)vinputs;
     2078
     2079        /* Get node coordinates and dof list: */
     2080        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     2081        GetDofList(&doflist[0],&numberofdofspernode);
     2082
     2083        /* Set Ke_gg to 0: */
     2084        for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;
     2085
     2086        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     2087        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     2088
     2089        /* Start  looping on the number of gaussian points: */
     2090        for (ig=0; ig<num_gauss; ig++){
     2091                /*Pick up the gaussian point: */
     2092                gauss_weight=*(gauss_weights+ig);
     2093                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     2094                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     2095                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     2096
     2097                /* Get Jacobian determinant: */
     2098                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     2099       
     2100                //Get L matrix if viscous basal drag present:
     2101                GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
     2102
     2103                DL_scalar=gauss_weight*Jdet;
     2104               
     2105                /*  Do the triple producte tL*D*L: */
     2106                TripleMultiply( &L[0][0],1,3,1,
     2107                                        &DL_scalar,1,1,0,
     2108                                        &L[0][0],1,3,0,
     2109                                        &Ke_gg_gaussian[0][0],0);
     2110
     2111               
     2112                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     2113                for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     2114        } //for (ig=0; ig<num_gauss; ig++
     2115
     2116        /*Add Ke_gg to global matrix Kgg: */
     2117        MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
     2118               
     2119        cleanup_and_return:
     2120        xfree((void**)&first_gauss_area_coord);
     2121        xfree((void**)&second_gauss_area_coord);
     2122        xfree((void**)&third_gauss_area_coord);
     2123        xfree((void**)&gauss_weights);
     2124}
     2125
    21042126#undef __FUNCT__
    21052127#define __FUNCT__ "Tria::CreateKMatrixDiagnosticSurfaceVert"
    2106 void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,ParameterInputs* inputs,int analysis_type){
     2128void  Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* vinputs,int analysis_type){
    21072129
    21082130        int i,j;
    … …  
    21442166        double Ke_gg[numdofs][numdofs]; //local element stiffness matrix
    21452167        double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix evaluated at the gaussian point.
     2168
     2169        ParameterInputs* inputs=NULL;
     2170
     2171        /*recover pointers: */
     2172        inputs=(ParameterInputs*)vinputs;
    21462173       
    21472174        /* Get node coordinates and dof list: */
    … …  
    22272254#undef __FUNCT__
    22282255#define __FUNCT__ "Tria::CreatePVectorDiagnosticBaseVert"
    2229 void  Tria::CreatePVectorDiagnosticBaseVert(Vec pg,ParameterInputs* inputs,int analysis_type){
     2256void  Tria::CreatePVectorDiagnosticBaseVert(Vec pg,void* vinputs,int analysis_type){
    22302257
    22312258        int             i,j;
    … …  
    22632290        /*input parameters for structural analysis (diagnostic): */
    22642291        double* velocity_param=NULL;
    2265         double  vx_list[numgrids];
    2266         double  vy_list[numgrids];
     2292        double  vx_list[numgrids]={0,0,0};
     2293        double  vy_list[numgrids]={0,0,0};
    22672294        double  vx,vy;
    22682295        double  meltingvalue;
    22692296        double  slope[2];
    22702297        double  dbdx,dbdy;
     2298        int     dofs1[1]={0};
     2299        int     dofs2[1]={1};
     2300
     2301        ParameterInputs* inputs=NULL;
     2302
     2303        /*recover pointers: */
     2304        inputs=(ParameterInputs*)vinputs;
    22712305
    22722306        /* recover input parameters: */
    2273         velocity_param=ParameterInputsRecover(inputs,"velocity");
    2274         if(!velocity_param)throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity");
     2307        if(!inputs->Recover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity");
     2308            inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes);
    22752309
    22762310        /* Get node coordinates and dof list: */
    … …  
    22812315        for(i=0;i<numdofs;i++) pe_g[i]=0.0;
    22822316
    2283         /*Initialize vx_list and vy_list: */
    2284         for(i=0;i<numgrids;i++){
    2285                 if(velocity_param){
    2286                         /*doflist is a list of dofs for the node. We have 1 dof per node for the vertical analysis. But the horizontal velocity
    2287                          *runs with 2 dofs per node! Therefore, we multiply by 2 the dof count, to extract the velocity from the inputs: */
    2288                         vx_list[i]=velocity_param[2*doflist[i*numberofdofspernode+0]+0];
    2289                         vy_list[i]=velocity_param[2*doflist[i*numberofdofspernode+0]+1];
    2290                 }
    2291                 else{
    2292                         vx_list[i]=0;
    2293                         vy_list[i]=0;
    2294                 }
    2295         }
    2296                
    22972317        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    22982318        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
Note: See TracChangeset for help on using the changeset viewer.