Ignore:
Timestamp:
04/01/10 14:10:22 (15 years ago)
Author:
Mathieu Morlighem
Message:

Fixed DG

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r3359 r3371  
    238238}
    239239/*}}}*/
     240/*FUNCTION Numericalflux::GetB {{{1*/
     241void Numericalflux::GetB(double* B, double gauss_coord){
     242
     243        const int numgrids=4;
     244        double l1l4[numgrids];
     245
     246        /*Get nodal functions: */
     247        GetNodalFunctions(&l1l4[0],gauss_coord);
     248
     249        /*Build B: */
     250        B[0] = +l1l4[0];
     251        B[1] = +l1l4[1];
     252        B[2] = -l1l4[0];
     253        B[3] = -l1l4[1];
     254}
     255/*}}}*/
    240256/*FUNCTION Numericalflux::GetDofList{{{1*/
    241257
    … …  
    254270                }
    255271        }
    256         else if (strcmp(type,"boundary")){
     272        else if (strcmp(type,"boundary")==0){
    257273                for(i=0;i<2;i++){
    258274                        nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
    … …  
    262278                }
    263279        }
    264         else ISSMERROR("type not supported yet");
     280        else ISSMERROR(exprintf("type %s not supported yet",type));
    265281
    266282        /*Assign output pointers:*/
    … …  
    288304}
    289305/*}}}*/
     306/*FUNCTION Numericalflux::GetL {{{1*/
     307void Numericalflux::GetL(double* L, double gauss_coord, int numdof){
     308
     309        const int numgrids=4;
     310        double l1l4[numgrids];
     311
     312        /*Check numdof*/
     313        ISSMASSERT(numdof==2 || numdof==4);
     314
     315        /*Get nodal functions: */
     316        GetNodalFunctions(&l1l4[0],gauss_coord);
     317
     318        /*Luild L: */
     319        L[0] = l1l4[0];
     320        L[1] = l1l4[1];
     321        if (numdof==4){
     322                L[2] = l1l4[2];
     323                L[3] = l1l4[3];
     324        }
     325}
     326/*}}}*/
    290327/*FUNCTION Numericalflux::GetName {{{1*/
    291328char* Numericalflux::GetName(void){
    … …  
    300337        l1l4[0]=-0.5*gauss_coord+0.5;
    301338        l1l4[1]=+0.5*gauss_coord+0.5;
    302         l1l4[0]=-0.5*gauss_coord+0.5;
    303         l1l4[1]=+0.5*gauss_coord+0.5;
    304 
     339        l1l4[2]=-0.5*gauss_coord+0.5;
     340        l1l4[3]=+0.5*gauss_coord+0.5;
     341
     342}
     343/*}}}*/
     344/*FUNCTION Numericalflux::GetNormal {{{1*/
     345void Numericalflux:: GetNormal(double* normal,double xyz_list[4][3]){
     346
     347        /*Build unit outward pointing vector*/
     348        const int numgrids=4;
     349        double vector[2];
     350        double norm;
     351
     352        vector[0]=xyz_list[1][0] - xyz_list[0][0];
     353        vector[1]=xyz_list[1][1] - xyz_list[0][1];
     354
     355        norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
     356
     357        normal[0]= + vector[1]/norm;
     358        normal[1]= - vector[0]/norm;
     359}
     360/*}}}*/
     361/*FUNCTION Numericalflux::GetParameterValue {{{1*/
     362void Numericalflux::GetParameterValue(double* pp, double* plist, double gauss_coord){
     363
     364        /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian
     365         * point specifie by gauss_l1l2l3: */
     366
     367        /*nodal functions: */
     368        double l1l4[4];
     369
     370        /*output: */
     371        double p;
     372
     373        GetNodalFunctions(&l1l4[0],gauss_coord);
     374
     375        p=l1l4[0]*plist[0]+l1l4[1]*plist[1];
     376
     377        /*Assign output pointers:*/
     378        *pp=p;
    305379}
    306380/*}}}*/
    … …  
    326400/*}}}*/
    327401/*FUNCTION Numericalflux::PenaltyCreateKMatrixInternal {{{1*/
    328 void  Numericalflux::PenaltyCreateKMatrixInternal(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
    329 
    330         ISSMERROR("not supported yet");
    331 
     402void  Numericalflux::PenaltyCreateKMatrixInternal(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
     403
     404        /* local declarations */
     405        int             i,j;
     406
     407        /* node data: */
     408        const int numgrids=4;
     409        const int NDOF1=1;
     410        const int numdof=NDOF1*numgrids;
     411        double    xyz_list[numgrids][3];
     412        double    normal[2];
     413        int       doflist[numdof];
     414        int       numberofdofspernode;
     415
     416        /* gaussian points: */
     417        int     num_gauss,ig;
     418        double* gauss_coords =NULL;
     419        double  gauss_coord;
     420        double* gauss_weights=NULL;
     421        double  gauss_weight;
     422
     423        /* matrices: */
     424        double B[numgrids];
     425        double L[numgrids];
     426        double DL1,DL2;
     427        double Ke_gg1[numdof][numdof];
     428        double Ke_gg2[numdof][numdof];
     429        double Ke_gg[numdof][numdof]={0.0};
     430        double Jdet;
     431
     432        /*input parameters for structural analysis (diagnostic): */
     433        double vx_list[numgrids]={0.0};
     434        double vy_list[numgrids]={0.0};
     435        double vx,vy;
     436        double UdotN;
     437        double dt;
     438        int    dofs[1]={0};
     439        int    found;
     440
     441        ParameterInputs* inputs=NULL;
     442
     443        /*recover pointers: */
     444        inputs=(ParameterInputs*)vinputs;
     445
     446        /*recover extra inputs from users, at current convergence iteration: */
     447        found=inputs->Recover("vx_average",&vx_list[0],1,dofs,numgrids,(void**)nodes);
     448        if(!found)ISSMERROR(" could not find vx_average in inputs!");
     449        found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes);
     450        if(!found)ISSMERROR(" could not find vy_average in inputs!");
     451        found=inputs->Recover("dt",&dt);
     452        if(!found)ISSMERROR(" could not find dt in inputs!");
     453
     454        /* Get node coordinates, dof list and normal vector: */
     455        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     456        GetDofList(&doflist[0],&numberofdofspernode);
     457        GetNormal(&normal[0],xyz_list);
     458
     459        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     460        num_gauss=2;
     461        GaussSegment(&gauss_coords, &gauss_weights,num_gauss);
     462
     463        /* Start  looping on the number of gaussian points: */
     464        for (ig=0; ig<num_gauss; ig++){
     465
     466                /*Pick up the gaussian point: */
     467                gauss_weight=gauss_weights[ig];
     468                gauss_coord =gauss_coords[ig];
     469
     470                /* Get Jacobian determinant: */
     471                GetJacobianDeterminant(&Jdet,xyz_list,gauss_coord);
     472
     473                //Get vx, vy and v.n
     474                GetParameterValue(&vx, &vx_list[0],gauss_coord);
     475                GetParameterValue(&vy, &vy_list[0],gauss_coord);
     476                UdotN=vx*normal[0]+vy*normal[1];
     477
     478                /*Get L and B: */
     479                GetL(&L[0],gauss_coord,numdof);
     480                GetB(&B[0],gauss_coord);
     481
     482                /*Compute DLs*/
     483                DL1=gauss_weight*Jdet*dt*UdotN/2;
     484                DL2=gauss_weight*Jdet*dt*fabs(UdotN)/2;
     485
     486                /*  Do the triple products: */
     487                TripleMultiply(&B[0],1,numdof,1,
     488                                        &DL1,1,1,0,
     489                                        &L[0],1,numdof,0,
     490                                        &Ke_gg1[0][0],0);
     491                TripleMultiply(&B[0],1,numdof,1,
     492                                        &DL2,1,1,0,
     493                                        &B[0],1,numdof,0,
     494                                        &Ke_gg2[0][0],0);
     495
     496                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     497                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg1[i][j];
     498                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg2[i][j];
     499
     500        } // for (ig=0; ig<num_gauss; ig++)
     501
     502        /*Add Ke_gg to global matrix Kgg: */
     503        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     504
     505
     506        xfree((void**)&gauss_coords);
     507        xfree((void**)&gauss_weights);
    332508}
    333509/*}}}*/
    334510/*FUNCTION Numericalflux::PenaltyCreateKMatrixBoundary {{{1*/
    335 void  Numericalflux::PenaltyCreateKMatrixBoundary(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
    336 
    337         ISSMERROR("type not supported yet");
     511void  Numericalflux::PenaltyCreateKMatrixBoundary(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
     512
     513        /* local declarations */
     514        int             i,j;
     515
     516        /* node data: */
     517        const int numgrids=2;
     518        const int NDOF1=1;
     519        const int numdof=NDOF1*numgrids;
     520        double    xyz_list[numgrids][3];
     521        double    normal[2];
     522        int       doflist[numdof];
     523        int       numberofdofspernode;
     524
     525        /* gaussian points: */
     526        int     num_gauss,ig;
     527        double* gauss_coords =NULL;
     528        double  gauss_coord;
     529        double* gauss_weights=NULL;
     530        double  gauss_weight;
     531
     532        /* matrices: */
     533        double L[numgrids];
     534        double DL;
     535        double Ke_gg[numdof][numdof]={0.0};
     536        double Ke_gg_gaussian[numdof][numdof];
     537        double Jdet;
     538
     539        /*input parameters for structural analysis (diagnostic): */
     540        double vx_list[numgrids]={0.0};
     541        double vy_list[numgrids]={0.0};
     542        double vx,vy;
     543        double mean_vx;
     544        double mean_vy;
     545        double UdotN;
     546        double dt;
     547        int    dofs[1]={0};
     548        int    found;
     549
     550        ParameterInputs* inputs=NULL;
     551
     552        /*recover pointers: */
     553        inputs=(ParameterInputs*)vinputs;
     554
     555        /*recover extra inputs from users, at current convergence iteration: */
     556        found=inputs->Recover("vx_average",&vx_list[0],1,dofs,numgrids,(void**)nodes);
     557        if(!found)ISSMERROR(" could not find vx_average in inputs!");
     558        found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes);
     559        if(!found)ISSMERROR(" could not find vy_average in inputs!");
     560        found=inputs->Recover("dt",&dt);
     561        if(!found)ISSMERROR(" could not find dt in inputs!");
     562
     563        /* Get node coordinates, dof list and normal vector: */
     564        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     565        GetDofList(&doflist[0],&numberofdofspernode);
     566        GetNormal(&normal[0],xyz_list);
     567
     568        /*Check wether it is an inflow or outflow BC*/
     569        mean_vx=(vx_list[0]+vx_list[1])/2.0;
     570        mean_vy=(vy_list[0]+vy_list[1])/2.0;
     571        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
     572        if (UdotN<=0){
     573                /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
     574                return;
     575        }
     576
     577        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     578        num_gauss=2;
     579        GaussSegment(&gauss_coords, &gauss_weights,num_gauss);
     580
     581        /* Start  looping on the number of gaussian points: */
     582        for (ig=0; ig<num_gauss; ig++){
     583
     584                /*Pick up the gaussian point: */
     585                gauss_weight=gauss_weights[ig];
     586                gauss_coord =gauss_coords[ig];
     587
     588                /* Get Jacobian determinant: */
     589                GetJacobianDeterminant(&Jdet,xyz_list,gauss_coord);
     590
     591                //Get vx, vy and v.n
     592                GetParameterValue(&vx, &vx_list[0],gauss_coord);
     593                GetParameterValue(&vy, &vy_list[0],gauss_coord);
     594                UdotN=vx*normal[0]+vy*normal[1];
     595
     596                /*Get L*/
     597                GetL(&L[0],gauss_coord,numdof);
     598
     599                /*Compute DLs*/
     600                DL=gauss_weight*Jdet*dt*UdotN;
     601
     602                /*Do triple product*/
     603                TripleMultiply(&L[0],1,numdof,1,
     604                                        &DL,1,1,0,
     605                                        &L[0],1,numdof,0,
     606                                        &Ke_gg_gaussian[0][0],0);
     607
     608                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     609                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     610
     611        } // for (ig=0; ig<num_gauss; ig++)
     612
     613        /*Add Ke_gg to global matrix Kgg: */
     614        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     615
     616        xfree((void**)&gauss_coords);
     617        xfree((void**)&gauss_weights);
    338618
    339619}
    … …  
    362642/*}}}*/
    363643/*FUNCTION Numericalflux::PenaltyCreatePVectorBoundary{{{1*/
    364 void  Numericalflux::PenaltyCreatePVectorBoundary(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){
    365 
    366         ISSMERROR("type not supported yet");
     644void  Numericalflux::PenaltyCreatePVectorBoundary(Vec pg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){
     645
     646        /* local declarations */
     647        int             i,j;
     648
     649        /* node data: */
     650        const int numgrids=2;
     651        const int NDOF1=1;
     652        const int numdof=NDOF1*numgrids;
     653        double    xyz_list[numgrids][3];
     654        double    normal[2];
     655        int       doflist[numdof];
     656        int       numberofdofspernode;
     657
     658        /* gaussian points: */
     659        int     num_gauss,ig;
     660        double* gauss_coords =NULL;
     661        double  gauss_coord;
     662        double* gauss_weights=NULL;
     663        double  gauss_weight;
     664
     665        /* matrices: */
     666        double L[numgrids];
     667        double DL;
     668        double pe_g[numdof]={0.0};
     669        double Jdet;
     670
     671        /*input parameters for structural analysis (diagnostic): */
     672        double vx_list[numgrids]={0.0};
     673        double vy_list[numgrids]={0.0};
     674        double vx,vy;
     675        double mean_vx;
     676        double mean_vy;
     677        double UdotN;
     678        double dt;
     679        int    dofs[1]={0};
     680        int    found;
     681
     682        ParameterInputs* inputs=NULL;
     683
     684        /*recover pointers: */
     685        inputs=(ParameterInputs*)vinputs;
     686
     687        /*recover extra inputs from users, at current convergence iteration: */
     688        found=inputs->Recover("vx_average",&vx_list[0],1,dofs,numgrids,(void**)nodes);
     689        if(!found)ISSMERROR(" could not find vx_average in inputs!");
     690        found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes);
     691        if(!found)ISSMERROR(" could not find vy_average in inputs!");
     692        found=inputs->Recover("dt",&dt);
     693        if(!found)ISSMERROR(" could not find dt in inputs!");
     694
     695        /* Get node coordinates, dof list and normal vector: */
     696        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     697        GetDofList(&doflist[0],&numberofdofspernode);
     698        GetNormal(&normal[0],xyz_list);
     699
     700        /*Check wether it is an inflow or outflow BC*/
     701        mean_vx=(vx_list[0]+vx_list[1])/2.0;
     702        mean_vy=(vy_list[0]+vy_list[1])/2.0;
     703        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
     704        if (UdotN>0){
     705                /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
     706                return;
     707        }
     708
     709        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     710        num_gauss=2;
     711        GaussSegment(&gauss_coords, &gauss_weights,num_gauss);
     712
     713        /* Start  looping on the number of gaussian points: */
     714        for (ig=0; ig<num_gauss; ig++){
     715
     716                /*Pick up the gaussian point: */
     717                gauss_weight=gauss_weights[ig];
     718                gauss_coord =gauss_coords[ig];
     719
     720                /* Get Jacobian determinant: */
     721                GetJacobianDeterminant(&Jdet,xyz_list,gauss_coord);
     722
     723                //Get vx, vy and v.n
     724                GetParameterValue(&vx, &vx_list[0],gauss_coord);
     725                GetParameterValue(&vy, &vy_list[0],gauss_coord);
     726                UdotN=vx*normal[0]+vy*normal[1];
     727
     728                /*Get L*/
     729                GetL(&L[0],gauss_coord,numdof);
     730
     731                /*Compute DL*/
     732                DL= - gauss_weight*Jdet*dt*UdotN;
     733
     734                /* Add value into pe_g: */
     735                for( i=0; i<numdof; i++) pe_g[i] += DL* L[i];
     736
     737        } // for (ig=0; ig<num_gauss; ig++)
     738
     739        /*Add pe_g to global matrix Kgg: */
     740        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
     741
     742        xfree((void**)&gauss_coords);
     743        xfree((void**)&gauss_weights);
    367744
    368745}
Note: See TracChangeset for help on using the changeset viewer.