Changeset 5741


Ignore:
Timestamp:
09/10/10 09:48:23 (15 years ago)
Author:
Mathieu Morlighem
Message:

removed some old gauss

Location:
issm/trunk/src/c/objects
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r5738 r5741  
    26052605
    26062606                }
    2607         } // for (ig=0; ig<num_gauss; ig++)
     2607        }
    26082608
    26092609        /*Add Ke_gg to global matrix Kgg: */
    … …  
    26922692                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg2[i][j];
    26932693
    2694         } // for (ig=0; ig<num_gauss; ig++)
     2694        }
    26952695
    26962696        /*Add Ke_gg to global matrix Kgg: */
    … …  
    28592859                }
    28602860
    2861         } // for (ig=0; ig<num_gauss; ig++)
     2861        }
    28622862
    28632863        /*Add Ke_gg to global matrix Kgg: */
    … …  
    29822982                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    29832983
    2984         } // for (ig=0; ig<num_gauss; ig++)
     2984        }
    29852985
    29862986        /*Add Ke_gg to global matrix Kgg: */
    … …  
    31093109                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    31103110
    3111         } // for (ig=0; ig<num_gauss; ig++)
     3111        }
    31123112
    31133113        /*Add Ke_gg to global matrix Kgg: */
    … …  
    32313231                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    32323232
    3233         } // for (ig=0; ig<num_gauss; ig++)
     3233        }
    32343234
    32353235        /*Add Ke_gg to global matrix Kgg: */
    … …  
    33513351                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    33523352
    3353         } // for (ig=0; ig<num_gauss; ig++)
     3353        }
    33543354
    33553355        /*Add Ke_gg to global matrix Kgg: */
    … …  
    34883488
    34893489
    3490         } //for (ig=0; ig<num_gauss; ig++)
     3490        }
    34913491
    34923492        /*Add Ke_gg to global matrix Kgg: */
    … …  
    37253725                }
    37263726
    3727         } // for (ig=0; ig<num_gauss; ig++)
     3727        }
    37283728
    37293729        /*Add Ke_gg to global matrix Kgg: */
    … …  
    38353835                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg2[i][j];
    38363836
    3837         } // for (ig=0; ig<num_gauss; ig++)
     3837        }
    38383838
    38393839        /*Add Ke_gg to global matrix Kgg: */
    … …  
    40344034                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g-dhdt_g)*L[i];
    40354035
    4036         } // for (ig=0; ig<num_gauss; ig++)
     4036        }
    40374037
    40384038        /*Add pe_g to global matrix Kgg: */
    … …  
    41044104                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g-dhdt_g)*L[i];
    41054105
    4106         } // for (ig=0; ig<num_gauss; ig++)
     4106        }
    41074107
    41084108        /*Add pe_g to global matrix Kgg: */
    … …  
    41704170                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g)*L[i];
    41714171
    4172         } // for (ig=0; ig<num_gauss; ig++)
     4172        }
    41734173
    41744174        /*Add pe_g to global matrix Kgg: */
    … …  
    42754275void Tria::CreatePVectorDiagnosticMacAyeal( Vec pg){
    42764276
    4277         int             i,j;
     4277        int             i,j,ig;
    42784278
    42794279        /* node data: */
    4280         const int    numdof=2*NUMVERTICES;
     4280        const int    numdof=NDOF2*NUMVERTICES;
    42814281        double       xyz_list[NUMVERTICES][3];
    42824282        int*         doflist=NULL;
    … …  
    42844284        /* parameters: */
    42854285        double  plastic_stress;
    4286         double  slope[NDOF2];
     4286        double  slope[2];
    42874287        double  driving_stress_baseline;
    4288 
    4289         /* gaussian points: */
    4290         int     num_gauss,ig;
    4291         double* first_gauss_area_coord  =  NULL;
    4292         double* second_gauss_area_coord =  NULL;
    4293         double* third_gauss_area_coord  =  NULL;
    4294         double* gauss_weights           =  NULL;
    4295         double  gauss_weight;
    4296         double  gauss_l1l2l3[3];
    4297 
    4298         /* Jacobian: */
     4288        GaussTria* gauss=NULL;
     4289
    42994290        double Jdet;
    4300 
    4301         /*nodal functions: */
    43024291        double l1l2l3[3];
    4303 
    4304         /*element vector at the gaussian points: */
    43054292        double  pe_g[numdof]={0.0};
    43064293        double  pe_g_gaussian[numdof];
    4307 
    4308         /*input parameters for structural analysis (diagnostic): */
    43094294        double  thickness;
    4310 
    4311         /*inputs: */
    43124295        bool onwater;
    43134296        int  drag_type;
    4314         Input* thickness_input=NULL;
    4315         Input* surface_input=NULL;
    4316         Input* drag_input=NULL;
    43174297
    43184298        /*retrieve inputs :*/
    43194299        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
    43204300        inputs->GetParameterValue(&drag_type,DragTypeEnum);
    4321         thickness_input=inputs->GetInput(ThicknessEnum);
    4322         surface_input=inputs->GetInput(SurfaceEnum);
    4323         drag_input=inputs->GetInput(DragCoefficientEnum);
     4301        Input* thickness_input=inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
     4302        Input* surface_input=inputs->GetInput(SurfaceEnum);     ISSMASSERT(surface_input);
     4303        Input* drag_input=inputs->GetInput(DragCoefficientEnum);ISSMASSERT(drag_input);
    43244304
    43254305        /*First, if we are on water, return empty vector: */
    4326         if(onwater)return;
     4306        if(onwater) return;
    43274307
    43284308        /* Get node coordinates and dof list: */
    … …  
    43304310        GetDofList(&doflist,MacAyealApproximationEnum);
    43314311
    4332 
    4333         /* Get gaussian points and weights: */
    4334         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/
    4335 
    43364312        /* Start  looping on the number of gaussian points: */
    4337         for (ig=0; ig<num_gauss; ig++){
    4338                 /*Pick up the gaussian point: */
    4339                 gauss_weight=*(gauss_weights+ig);
    4340                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    4341                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    4342                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    4343 
    4344                 /*Compute thickness at gaussian point: */
    4345                 thickness_input->GetParameterValue(&thickness, &gauss_l1l2l3[0]);
    4346                 surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
     4313        gauss=new GaussTria(2);
     4314        for(ig=gauss->begin();ig<gauss->end();ig++){
     4315
     4316                gauss->GaussPoint(ig);
     4317
     4318                thickness_input->GetParameterValue(&thickness,gauss);
     4319                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    43474320               
    43484321                /*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the
    43494322                 * element itself: */
    43504323                if(drag_type==1){
    4351                         drag_input->GetParameterValue(&plastic_stress, &gauss_l1l2l3[0]);
     4324                        drag_input->GetParameterValue(&plastic_stress,gauss);
    43524325                }
    43534326
    43544327                /* Get Jacobian determinant: */
    4355                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     4328                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    43564329               
    43574330                 /*Get nodal functions: */
    4358                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     4331                GetNodalFunctions(l1l2l3, gauss);
    43594332
    43604333                /*Compute driving stress: */
    … …  
    43654338                        for (i=0;i<NUMVERTICES;i++){
    43664339                                for (j=0;j<NDOF2;j++){
    4367                                         pe_g_gaussian[i*NDOF2+j]=(-driving_stress_baseline*slope[j]-plastic_stress)*Jdet*gauss_weight*l1l2l3[i];
     4340                                        pe_g_gaussian[i*NDOF2+j]=(-driving_stress_baseline*slope[j]-plastic_stress)*Jdet*gauss->weight*l1l2l3[i];
    43684341                                }
    43694342                        }
    … …  
    43724345                        for (i=0;i<NUMVERTICES;i++){
    43734346                                for (j=0;j<NDOF2;j++){
    4374                                         pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss_weight*l1l2l3[i];
     4347                                        pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*l1l2l3[i];
    43754348                                }
    43764349                        }
    … …  
    43804353                for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
    43814354
    4382         } //for (ig=0; ig<num_gauss; ig++)
     4355        }
    43834356
    43844357        /*Add pe_g to global vector pg: */
    … …  
    43864359
    43874360        /*Free ressources:*/
    4388         xfree((void**)&first_gauss_area_coord);
    4389         xfree((void**)&second_gauss_area_coord);
    4390         xfree((void**)&third_gauss_area_coord);
    4391         xfree((void**)&gauss_weights);
     4361        delete gauss;
    43924362        xfree((void**)&doflist);
    43934363}
    … …  
    50435013                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i];
    50445014
    5045         } // for (ig=0; ig<num_gauss; ig++)
     5015        }
    50465016
    50475017        /*Add pe_g to global matrix Kgg: */
    … …  
    51165086                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i];
    51175087
    5118         } // for (ig=0; ig<num_gauss; ig++)
     5088        }
    51195089
    51205090        /*Add pe_g to global matrix Kgg: */
    … …  
    51955165                for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
    51965166
    5197         } //for (ig=0; ig<num_gauss; ig++)
     5167        }
    51985168
    51995169        /*Add pe_g to global vector pg: */
    … …  
    57665736void  Tria::GradjDragStokes(Vec gradient){
    57675737
    5768         int i;
     5738        int i,ig;
    57695739
    57705740        /* node data: */
    … …  
    57775747        double alpha_complement;
    57785748        Friction* friction=NULL;
    5779 
    5780         /* gaussian points: */
    5781         int     num_gauss,ig;
    5782         double* first_gauss_area_coord  =  NULL;
    5783         double* second_gauss_area_coord =  NULL;
    5784         double* third_gauss_area_coord  =  NULL;
    5785         double* gauss_weights           =  NULL;
    5786         double  gauss_weight;
    5787         double  gauss_l1l2l3[3];
     5749        GaussTria* gauss=NULL;
    57885750
    57895751        /* parameters: */
    … …  
    57985760        double  grade_g[NUMVERTICES]={0.0};
    57995761        double  grade_g_gaussian[NUMVERTICES];
    5800 
    5801         /* Jacobian: */
    58025762        double Jdet;
    5803 
    5804         /*nodal functions: */
    58055763        double l1l2l3[3];
    5806 
    5807         /* strain rate: */
    58085764        double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    5809 
    5810         /*inputs: */
    58115765        bool shelf;
    58125766        int  drag_type;
    5813 
    5814         /*parameters: */
    58155767        double  cm_noisedmp;
    5816 
    58175768        int analysis_type;
    58185769
    … …  
    58235774        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
    58245775        inputs->GetParameterValue(&drag_type,DragTypeEnum);
     5776        Input* drag_input=inputs->GetInput(DragCoefficientEnum); ISSMASSERT(drag_input);
    58255777
    58265778        /*retrieve some parameters: */
    … …  
    58385790        friction=new Friction("2d",inputs,matpar,analysis_type);
    58395791
    5840         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    5841         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
    5842 
    58435792        /* Start  looping on the number of gaussian points: */
    5844         for (ig=0; ig<num_gauss; ig++){
    5845                 /*Pick up the gaussian point: */
    5846                 gauss_weight=*(gauss_weights+ig);
    5847                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    5848                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    5849                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     5793        gauss=new GaussTria(4);
     5794        for(ig=gauss->begin();ig<gauss->end();ig++){
     5795
     5796                gauss->GaussPoint(ig);
    58505797
    58515798                /*Recover alpha_complement and drag: */
    5852                 if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxEnum,VyEnum);
     5799                if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum);
    58535800                else alpha_complement=0;
    5854                 inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum);
     5801                inputs->GetParameterValue(&drag,gauss,DragCoefficientEnum);
    58555802
    58565803                /*recover lambda mu and xi: */
    5857                 inputs->GetParameterValue(&lambda, &gauss_l1l2l3[0],AdjointxEnum);
    5858                 inputs->GetParameterValue(&mu, &gauss_l1l2l3[0],AdjointyEnum);
    5859                 inputs->GetParameterValue(&xi, &gauss_l1l2l3[0],AdjointzEnum);
     5804                inputs->GetParameterValue(&lambda,gauss,AdjointxEnum);
     5805                inputs->GetParameterValue(&mu, gauss,AdjointyEnum);
     5806                inputs->GetParameterValue(&xi, gauss,AdjointzEnum);
    58605807
    58615808                /*recover vx vy and vz: */
    5862                 inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum);
    5863                 inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum);
    5864                 inputs->GetParameterValue(&vz, &gauss_l1l2l3[0],VzEnum);
     5809                inputs->GetParameterValue(&vx, gauss,VxEnum);
     5810                inputs->GetParameterValue(&vy, gauss,VyEnum);
     5811                inputs->GetParameterValue(&vz, gauss,VzEnum);
    58655812
    58665813                /*Get normal vector to the bed */
    … …  
    58725819
    58735820                /* Get Jacobian determinant: */
    5874                 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     5821                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss);
    58755822
    58765823                /* Get nodal functions value at gaussian point:*/
    5877                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     5824                GetNodalFunctions(l1l2l3, gauss);
    58785825
    58795826                /*Get nodal functions derivatives*/
    5880                 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
     5827                GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
    58815828
    58825829                /*Get k derivative: dk/dx */
    5883                 inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);
     5830                drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
    58845831
    58855832                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    … …  
    58905837                                                -mu    *(2*drag*alpha_complement*(vy - vz*bed_normal[1]*bed_normal[2]))
    58915838                                                -xi    *(2*drag*alpha_complement*(-vx*bed_normal[0]*bed_normal[2]-vy*bed_normal[1]*bed_normal[2]))
    5892                                                 )*Jdet*gauss_weight*l1l2l3[i];
     5839                                                )*Jdet*gauss->weight*l1l2l3[i];
    58935840
    58945841                        //Add regularization term
    5895                         grade_g_gaussian[i]+= - cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
     5842                        grade_g_gaussian[i]+= - cm_noisedmp*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
    58965843                }
    58975844
    … …  
    59035850        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
    59045851
    5905         /*Add grade_g to the inputs of this element: */
    5906         this->inputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0]));
    5907 
    5908         xfree((void**)&first_gauss_area_coord);
    5909         xfree((void**)&second_gauss_area_coord);
    5910         xfree((void**)&third_gauss_area_coord);
    5911         xfree((void**)&gauss_weights);
    59125852        delete friction;
     5853        delete gauss;
    59135854
    59145855}

  • issm/trunk/src/c/objects/Loads/Numericalflux.cpp

    r5739 r5741  
    652652        int count=0;
    653653        int type;
    654         int numberofnodes=2;
    655 
    656         /*output: */
     654        int numberofnodes;
    657655        int* doflist=NULL;
    658656       
    659         /*recover type: */
    660         inputs->GetParameterValue(&type,TypeEnum);
    661 
    662657        /*Some checks for debugging*/
    663658        ISSMASSERT(nodes);
    664659               
    665660        /*How many nodes? :*/
    666         if(type==InternalEnum)numberofnodes=4;
    667         else if(type==BoundaryEnum)numberofnodes=2;
    668         else ISSMERROR("type %s not supported yet",type);
     661        inputs->GetParameterValue(&type,TypeEnum);
     662        switch(type){
     663                case InternalEnum:
     664                        numberofnodes=NUMVERTICES_INTERNAL; break;
     665                case BoundaryEnum:
     666                        numberofnodes=NUMVERTICES_BOUNDARY; break;
     667                default:
     668                        ISSMERROR("type %s not supported yet",type);
     669        }
    669670       
    670671        /*Figure out size of doflist: */
    … …  
    691692
    692693        /*Build unit outward pointing vector*/
    693         const int numgrids=4;
    694694        double vector[2];
    695695        double norm;
Note: See TracChangeset for help on using the changeset viewer.