Changeset 5639


Ignore:
Timestamp:
08/31/10 16:00:44 (15 years ago)
Author:
Mathieu Morlighem
Message:

Final touch of GaussTria

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r5638 r5639  
    46774677
    46784678        /* gaussian points: */
    4679         int     num_gauss,ig;
    4680         double* first_gauss_area_coord  =  NULL;
    4681         double* second_gauss_area_coord =  NULL;
    4682         double* third_gauss_area_coord  =  NULL;
    4683         double* gauss_weights           =  NULL;
    4684         double  gauss_weight;
    4685         double  gauss_l1l2l3[3];
     4679        int     ig;
     4680        GaussTria* gauss=NULL;
    46864681
    46874682        /* parameters: */
    … …  
    47124707        weights_input     =inputs->GetInput(WeightsEnum);ISSMASSERT(weights_input);
    47134708
    4714         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4715         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    4716 
    47174709        /* Start  looping on the number of gaussian points: */
    4718         for (ig=0; ig<num_gauss; ig++){
    4719                 /*Pick up the gaussian point: */
    4720                 gauss_weight=*(gauss_weights+ig);
    4721                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    4722                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    4723                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     4710        gauss=new GaussTria(2);
     4711        for(ig=gauss->begin();ig<gauss->end();ig++){
     4712
     4713                gauss->GaussPoint(ig);
    47244714
    47254715                /* Get Jacobian determinant: */
    4726                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     4716                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    47274717
    47284718                /* Get nodal functions value at gaussian point:*/
    4729                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     4719                GetNodalFunctions(l1l2l3, gauss);
    47304720
    47314721                /*Get parameters at gauss point*/
    4732                 thickness_input->GetParameterValue(&thickness, &gauss_l1l2l3[0]);
    4733                 thicknessobs_input->GetParameterValue(&thicknessobs, &gauss_l1l2l3[0]);
    4734                 weights_input->GetParameterValue(&weight, &gauss_l1l2l3[0]);
     4722                thickness_input->GetParameterValue(&thickness, gauss);
     4723                thicknessobs_input->GetParameterValue(&thicknessobs, gauss);
     4724                weights_input->GetParameterValue(&weight, gauss);
    47354725
    47364726                for (i=0;i<numvertices;i++){
    4737                         pe_g_gaussian[i]=(thicknessobs-thickness)*weight*Jdet*gauss_weight*l1l2l3[i];
     4727                        pe_g_gaussian[i]=(thicknessobs-thickness)*weight*Jdet*gauss->weight*l1l2l3[i];
    47384728                }
    47394729
    … …  
    47474737        VecSetValues(p_g,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    47484738
    4749         /*Free ressources:*/
    4750         xfree((void**)&first_gauss_area_coord);
    4751         xfree((void**)&second_gauss_area_coord);
    4752         xfree((void**)&third_gauss_area_coord);
    4753         xfree((void**)&gauss_weights);
     4739        /*Clean up and return*/
     4740        delete gauss;
    47544741        xfree((void**)&doflist);
    47554742}
    … …  
    47774764
    47784765        /* gaussian points: */
    4779         int     num_gauss,ig;
    4780         double* first_gauss_area_coord  =  NULL;
    4781         double* second_gauss_area_coord =  NULL;
    4782         double* third_gauss_area_coord  =  NULL;
    4783         double* gauss_weights           =  NULL;
    4784         double  gauss_weight;
    4785         double  gauss_l1l2l3[3];
     4766        int     ig;
     4767        GaussTria* gauss=NULL;
    47864768
    47874769        /* parameters: */
    … …  
    49454927        }
    49464928
    4947         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4948         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    4949 
    49504929        /* Start  looping on the number of gaussian points: */
    4951         for (ig=0; ig<num_gauss; ig++){
    4952                 /*Pick up the gaussian point: */
    4953                 gauss_weight=*(gauss_weights+ig);
    4954                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    4955                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    4956                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     4930        gauss=new GaussTria(2);
     4931        for(ig=gauss->begin();ig<gauss->end();ig++){
     4932
     4933                gauss->GaussPoint(ig);
    49574934
    49584935                /* Get Jacobian determinant: */
    4959                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     4936                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    49604937
    49614938                /* Get nodal functions value at gaussian point:*/
    4962                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     4939                GetNodalFunctions(l1l2l3, gauss);
    49634940
    49644941                /*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */
    49654942
    49664943                /*Compute absolute(x/y) at gaussian point: */
    4967                 TriaRef::GetParameterValue(&dux, &dux_list[0],gauss_l1l2l3);
    4968                 TriaRef::GetParameterValue(&duy, &duy_list[0],gauss_l1l2l3);
     4944                TriaRef::GetParameterValue(&dux, &dux_list[0],gauss);
     4945                TriaRef::GetParameterValue(&duy, &duy_list[0],gauss);
    49694946
    49704947                /*compute Du*/
    49714948                for (i=0;i<numvertices;i++){
    4972                         pe_g_gaussian[i*NDOF2+0]=dux*Jdet*gauss_weight*l1l2l3[i];
    4973                         pe_g_gaussian[i*NDOF2+1]=duy*Jdet*gauss_weight*l1l2l3[i];
     4949                        pe_g_gaussian[i*NDOF2+0]=dux*Jdet*gauss->weight*l1l2l3[i];
     4950                        pe_g_gaussian[i*NDOF2+1]=duy*Jdet*gauss->weight*l1l2l3[i];
    49744951                }
    49754952
    … …  
    49834960        VecSetValues(p_g,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    49844961
    4985         /*Free ressources:*/
    4986         xfree((void**)&first_gauss_area_coord);
    4987         xfree((void**)&second_gauss_area_coord);
    4988         xfree((void**)&third_gauss_area_coord);
    4989         xfree((void**)&gauss_weights);
     4962        /*Clean up and return*/
     4963        delete gauss;
    49904964        xfree((void**)&doflist);
    49914965}
    … …  
    50134987
    50144988        /* gaussian points: */
    5015         int     num_gauss,ig;
    5016         double* first_gauss_area_coord  =  NULL;
    5017         double* second_gauss_area_coord =  NULL;
    5018         double* third_gauss_area_coord  =  NULL;
    5019         double* gauss_weights           =  NULL;
    5020         double  gauss_weight;
    5021         double  gauss_l1l2l3[3];
     4989        int     ig;
     4990        GaussTria* gauss=NULL;
    50224991
    50234992        /* parameters: */
    … …  
    51815150        }
    51825151
    5183         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    5184         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    5185 
    51865152        /* Start  looping on the number of gaussian points: */
    5187         for (ig=0; ig<num_gauss; ig++){
    5188                 /*Pick up the gaussian point: */
    5189                 gauss_weight=*(gauss_weights+ig);
    5190                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    5191                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    5192                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     5153        gauss=new GaussTria(2);
     5154        for(ig=gauss->begin();ig<gauss->end();ig++){
     5155
     5156                gauss->GaussPoint(ig);
    51935157
    51945158                /* Get Jacobian determinant: */
    5195                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     5159                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    51965160
    51975161                /* Get nodal functions value at gaussian point:*/
    5198                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     5162                GetNodalFunctions(l1l2l3, gauss);
    51995163
    52005164                /*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */
    52015165
    52025166                /*Compute absolute(x/y) at gaussian point: */
    5203                 TriaRef::GetParameterValue(&dux, &dux_list[0],gauss_l1l2l3);
    5204                 TriaRef::GetParameterValue(&duy, &duy_list[0],gauss_l1l2l3);
     5167                TriaRef::GetParameterValue(&dux, &dux_list[0],gauss);
     5168                TriaRef::GetParameterValue(&duy, &duy_list[0],gauss);
    52055169
    52065170                /*compute Du*/
    52075171                for (i=0;i<numvertices;i++){
    5208                         pe_g[i*NDOF4+0]+=dux*Jdet*gauss_weight*l1l2l3[i];
    5209                         pe_g[i*NDOF4+1]+=duy*Jdet*gauss_weight*l1l2l3[i];
     5172                        pe_g[i*NDOF4+0]+=dux*Jdet*gauss->weight*l1l2l3[i];
     5173                        pe_g[i*NDOF4+1]+=duy*Jdet*gauss->weight*l1l2l3[i];
    52105174                }
    52115175        }
    … …  
    52145178        VecSetValues(p_g,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    52155179
    5216         /*Free ressources:*/
    5217         xfree((void**)&first_gauss_area_coord);
    5218         xfree((void**)&second_gauss_area_coord);
    5219         xfree((void**)&third_gauss_area_coord);
    5220         xfree((void**)&gauss_weights);
     5180        /*Clean up and return*/
     5181        delete gauss;
    52215182        xfree((void**)&doflist);
    52225183}
    … …  
    52445205        Input* slopey_input=NULL;
    52455206        Input* thickness_input=NULL;
    5246         double gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
     5207        GaussTria* gauss=NULL;
    52475208
    52485209        /*recover some inputs: */
    … …  
    52665227
    52675228        /*Spawn 3 sing elements: */
    5268         for(i=0;i<3;i++){
     5229        gauss=new GaussTria();
     5230        for(i=0;i<numvertices;i++){
     5231
     5232                gauss->GaussVertex(i);
     5233
    52695234                connectivity=nodes[i]->GetConnectivity();
    52705235
    5271                 thickness_input->GetParameterValue(&thickness,&gauss[i][0]);
    5272                 slopex_input->GetParameterValue(&slope[0],&gauss[i][0]);
    5273                 slopey_input->GetParameterValue(&slope[1],&gauss[i][0]);
     5236                thickness_input->GetParameterValue(&thickness,gauss);
     5237                slopex_input->GetParameterValue(&slope[0],gauss);
     5238                slopey_input->GetParameterValue(&slope[1],gauss);
    52745239                slope2=pow(slope[0],2)+pow(slope[1],2);
    52755240
    … …  
    52875252
    52885253        VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);
    5289        
    5290         /*Free ressources:*/
     5254
     5255        /*Clean up and return*/
     5256        delete gauss;
    52915257        xfree((void**)&doflist);
    52925258}
    … …  
    52945260/*FUNCTION Tria::CreatePVectorPrognostic_CG {{{1*/
    52955261void  Tria::CreatePVectorPrognostic_CG(Vec pg){
    5296 
    52975262
    52985263        /* local declarations */
    … …  
    53085273
    53095274        /* gaussian points: */
    5310         int     num_gauss,ig;
    5311         double* first_gauss_area_coord  =  NULL;
    5312         double* second_gauss_area_coord =  NULL;
    5313         double* third_gauss_area_coord  =  NULL;
    5314         double* gauss_weights           =  NULL;
    5315         double  gauss_weight;
    5316         double  gauss_l1l2l3[3];
     5275        int     ig;
     5276        GaussTria* gauss=NULL;
    53175277
    53185278        /* matrix */
    … …  
    53395299        GetDofList(&doflist);
    53405300
    5341         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    5342         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    5343 
    53445301        /*Retrieve all inputs we will be needing: */
    53455302        accumulation_input=inputs->GetInput(AccumulationRateEnum);
    … …  
    53485305
    53495306        /* Start  looping on the number of gaussian points: */
    5350         for (ig=0; ig<num_gauss; ig++){
    5351                 /*Pick up the gaussian point: */
    5352                 gauss_weight=*(gauss_weights+ig);
    5353                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    5354                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    5355                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     5307        gauss=new GaussTria(2);
     5308        for(ig=gauss->begin();ig<gauss->end();ig++){
     5309
     5310                gauss->GaussPoint(ig);
    53565311
    53575312                /* Get Jacobian determinant: */
    5358                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     5313                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    53595314
    53605315                /*Get L matrix: */
    5361                 GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     5316                GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode);
    53625317
    53635318                /* Get accumulation, melting and thickness at gauss point */
    5364                 accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]);
    5365                 melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]);
    5366                 thickness_input->GetParameterValue(&thickness_g, &gauss_l1l2l3[0]);
     5319                accumulation_input->GetParameterValue(&accumulation_g,gauss);
     5320                melting_input->GetParameterValue(&melting_g,gauss);
     5321                thickness_input->GetParameterValue(&thickness_g,gauss);
    53675322
    53685323                /* Add value into pe_g: */
    5369                 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i];
     5324                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i];
    53705325
    53715326        } // for (ig=0; ig<num_gauss; ig++)
    … …  
    53745329        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    53755330
    5376         /*Free ressources:*/
    5377         xfree((void**)&first_gauss_area_coord);
    5378         xfree((void**)&second_gauss_area_coord);
    5379         xfree((void**)&third_gauss_area_coord);
    5380         xfree((void**)&gauss_weights);
     5331        /*Clean up and return*/
     5332        delete gauss;
    53815333        xfree((void**)&doflist);
    53825334
    … …  
    53985350
    53995351        /* gaussian points: */
    5400         int     num_gauss,ig;
    5401         double* first_gauss_area_coord  =  NULL;
    5402         double* second_gauss_area_coord =  NULL;
    5403         double* third_gauss_area_coord  =  NULL;
    5404         double* gauss_weights           =  NULL;
    5405         double  gauss_weight;
    5406         double  gauss_l1l2l3[3];
     5352        int     ig;
     5353        GaussTria* gauss=NULL;
    54075354
    54085355        /* matrix */
    … …  
    54275374        GetDofList(&doflist);
    54285375
    5429         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    5430         GaussLegendreTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    5431 
    54325376        /*Retrieve all inputs we will be needing: */
    54335377        accumulation_input=inputs->GetInput(AccumulationRateEnum);
    … …  
    54365380
    54375381        /* Start  looping on the number of gaussian points: */
    5438         for (ig=0; ig<num_gauss; ig++){
    5439                 /*Pick up the gaussian point: */
    5440                 gauss_weight=*(gauss_weights+ig);
    5441                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    5442                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    5443                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     5382        gauss=new GaussTria(2);
     5383        for(ig=gauss->begin();ig<gauss->end();ig++){
     5384
     5385                gauss->GaussPoint(ig);
    54445386
    54455387                /* Get Jacobian determinant: */
    5446                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     5388                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    54475389
    54485390                /*Get L matrix: */
    5449                 GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     5391                GetL(&L[0], &xyz_list[0][0], gauss,numberofdofspernode);
    54505392
    54515393                /* Get accumulation, melting and thickness at gauss point */
    5452                 accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]);
    5453                 melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]);
    5454                 thickness_input->GetParameterValue(&thickness_g, &gauss_l1l2l3[0]);
     5394                accumulation_input->GetParameterValue(&accumulation_g,gauss);
     5395                melting_input->GetParameterValue(&melting_g,gauss);
     5396                thickness_input->GetParameterValue(&thickness_g,gauss);
    54555397
    54565398                /* Add value into pe_g: */
    5457                 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i];
     5399                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i];
    54585400
    54595401        } // for (ig=0; ig<num_gauss; ig++)
    … …  
    54625404        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    54635405
    5464         /*Free ressources:*/
    5465         xfree((void**)&first_gauss_area_coord);
    5466         xfree((void**)&second_gauss_area_coord);
    5467         xfree((void**)&third_gauss_area_coord);
    5468         xfree((void**)&gauss_weights);
     5406        /*Clean up and return*/
     5407        delete gauss;
    54695408        xfree((void**)&doflist);
    54705409}
    … …  
    54835422       
    54845423        /* gaussian points: */
    5485         int     num_gauss,ig;
    5486         double* first_gauss_area_coord  =  NULL;
    5487         double* second_gauss_area_coord =  NULL;
    5488         double* third_gauss_area_coord  =  NULL;
    5489         double* gauss_weights           =  NULL;
    5490         double  gauss_weight;
    5491         double  gauss_l1l2l3[3];
     5424        int     ig;
     5425        GaussTria* gauss=NULL;
    54925426
    54935427        /* Jacobian: */
    … …  
    55135447        GetDofList(&doflist);
    55145448
    5515         /* Get gaussian points and weights: */
    5516         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*/
    5517 
    55185449        /*Retrieve all inputs we will be needing: */
    55195450        if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==SurfaceSlopeYAnalysisEnum)){
    … …  
    55255456               
    55265457        /* Start  looping on the number of gaussian points: */
    5527         for (ig=0; ig<num_gauss; ig++){
    5528                 /*Pick up the gaussian point: */
    5529                 gauss_weight=*(gauss_weights+ig);
    5530                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    5531                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    5532                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    5533 
    5534                 slope_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]);
     5458        gauss=new GaussTria(2);
     5459        for(ig=gauss->begin();ig<gauss->end();ig++){
     5460
     5461                gauss->GaussPoint(ig);
     5462
     5463                slope_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    55355464
    55365465                /* Get Jacobian determinant: */
    5537                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     5466                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    55385467               
    55395468                 /*Get nodal functions: */
    5540                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     5469                GetNodalFunctions(l1l2l3, gauss);
    55415470
    55425471                /*Build pe_g_gaussian vector: */
    55435472                if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==BedSlopeXAnalysisEnum)){
    5544                         for(i=0;i<numdof;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[0]*l1l2l3[i];
     5473                        for(i=0;i<numdof;i++) pe_g_gaussian[i]=Jdet*gauss->weight*slope[0]*l1l2l3[i];
    55455474                }
    55465475                if ( (analysis_type==SurfaceSlopeYAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){
    5547                         for(i=0;i<numdof;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[1]*l1l2l3[i];
     5476                        for(i=0;i<numdof;i++) pe_g_gaussian[i]=Jdet*gauss->weight*slope[1]*l1l2l3[i];
    55485477                }
    55495478
    … …  
    55565485        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    55575486
    5558         /*Free ressources:*/
    5559         xfree((void**)&first_gauss_area_coord);
    5560         xfree((void**)&second_gauss_area_coord);
    5561         xfree((void**)&third_gauss_area_coord);
    5562         xfree((void**)&gauss_weights);
     5487        /*Clean up and return*/
     5488        delete gauss;
    55635489        xfree((void**)&doflist);
    55645490}
    … …  
    55885514
    55895515        /* gaussian points: */
    5590         int     num_area_gauss,ig;
    5591         double* gauss_weights  =  NULL;
    5592         double* first_gauss_area_coord  =  NULL;
    5593         double* second_gauss_area_coord =  NULL;
    5594         double* third_gauss_area_coord  =  NULL;
    5595         double  gauss_weight;
    5596         double  gauss_coord[3];
     5516        int     ig;
     5517        GaussTria* gauss=NULL;
    55975518
    55985519        /*matrices: */
    … …  
    56255546        /* Ice/ocean heat exchange flux on ice shelf base */
    56265547
    5627         GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    5628 
    56295548        /*Retrieve all inputs we will be needing: */
    56305549        pressure_input=inputs->GetInput(PressureEnum);
    56315550
    56325551        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    5633         for (ig=0; ig<num_area_gauss; ig++){
    5634                 gauss_weight=*(gauss_weights+ig);
    5635                 gauss_coord[0]=*(first_gauss_area_coord+ig);
    5636                 gauss_coord[1]=*(second_gauss_area_coord+ig);
    5637                 gauss_coord[2]=*(third_gauss_area_coord+ig);
     5552        gauss=new GaussTria(2);
     5553        for(ig=gauss->begin();ig<gauss->end();ig++){
     5554
     5555                gauss->GaussPoint(ig);
    56385556
    56395557                //Get the Jacobian determinant
    5640                 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss_coord);
     5558                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss);
    56415559
    56425560                /*Get nodal functions values: */
    5643                 GetNodalFunctions(&l1l2l3[0], gauss_coord);
     5561                GetNodalFunctions(&l1l2l3[0], gauss);
    56445562
    56455563                /*Get geothermal flux and basal friction */
    5646                 pressure_input->GetParameterValue(&pressure, &gauss_coord[0]);
     5564                pressure_input->GetParameterValue(&pressure,gauss);
    56475565                t_pmp=meltingpoint-beta*pressure;
    56485566
    56495567                /*Calculate scalar parameter*/
    5650                 scalar_ocean=gauss_weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice);
     5568                scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice);
    56515569                if(dt){
    56525570                        scalar_ocean=dt*scalar_ocean;
    … …  
    56615579        VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
    56625580
    5663         /*Free ressources:*/
    5664         xfree((void**)&first_gauss_area_coord);
    5665         xfree((void**)&second_gauss_area_coord);
    5666         xfree((void**)&third_gauss_area_coord);
    5667         xfree((void**)&gauss_weights);
     5581        /*Clean up and return*/
     5582        delete gauss;
    56685583        xfree((void**)&doflist);
    56695584}
    … …  
    56935608        double geothermalflux_value;
    56945609        double alpha2_list[numvertices];                                    //TO BE DELETED
    5695         double gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
     5610        double gaussgrids[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
    56965611        double vx_list[numvertices]; //TO BE DELETED
    56975612        double vy_list[numvertices]; //TO BE DELETED
    … …  
    56995614
    57005615        /* gaussian points: */
    5701         int     num_area_gauss,ig;
    5702         double* gauss_weights  =  NULL;
    5703         double* first_gauss_area_coord  =  NULL;
    5704         double* second_gauss_area_coord =  NULL;
    5705         double* third_gauss_area_coord  =  NULL;
    5706         double  gauss_weight;
    5707         double  gauss_coord[3];
     5616        int     ig;
     5617        GaussTria* gauss=NULL;
    57085618
    57095619        /*matrices: */
    … …  
    57485658        /*COMPUT alpha2_list and basalfriction_list (TO BE DELETED)*/
    57495659        for(i=0;i<numvertices;i++){
    5750                 friction->GetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum); //TO BE DELETED
    5751         }
    5752         vx_input->GetParameterValues(&vx_list[0],&gauss[0][0],3); //TO BE DELETED
    5753         vy_input->GetParameterValues(&vy_list[0],&gauss[0][0],3); //TO BE DELETED
     5660                friction->GetAlpha2(&alpha2_list[i],&gaussgrids[i][0],VxEnum,VyEnum,VzEnum); //TO BE DELETED
     5661        }
     5662        vx_input->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3); //TO BE DELETED
     5663        vy_input->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3); //TO BE DELETED
    57545664        for(i=0;i<numvertices;i++){
    57555665                basalfriction_list[i]=alpha2_list[i]*(pow(vx_list[i],(double)2.0)+pow(vy_list[i],(double)2.0)); //TO BE DELETED
    … …  
    57575667       
    57585668        /* Ice/ocean heat exchange flux on ice shelf base */
    5759         GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    57605669
    57615670        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    5762         for (ig=0; ig<num_area_gauss; ig++){
    5763                 gauss_weight=*(gauss_weights+ig);
    5764                 gauss_coord[0]=*(first_gauss_area_coord+ig);
    5765                 gauss_coord[1]=*(second_gauss_area_coord+ig);
    5766                 gauss_coord[2]=*(third_gauss_area_coord+ig);
     5671        gauss=new GaussTria(2);
     5672        for(ig=gauss->begin();ig<gauss->end();ig++){
     5673
     5674                gauss->GaussPoint(ig);
    57675675
    57685676                //Get the Jacobian determinant
    5769                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss_coord);
     5677                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss);
    57705678
    57715679                /*Get nodal functions values: */
    5772                 GetNodalFunctions(&l1l2l3[0], gauss_coord);
     5680                GetNodalFunctions(&l1l2l3[0], gauss);
    57735681
    57745682                /*Get geothermal flux and basal friction */
    5775                 geothermalflux_input->GetParameterValue(&geothermalflux_value, &gauss_coord[0]);
     5683                geothermalflux_input->GetParameterValue(&geothermalflux_value,gauss);
    57765684       
    57775685                /*Friction: */
    57785686                //friction->GetAlpha2(&alpha2,&gauss_coord[0],VxEnum,VyEnum,VzEnum);
    5779                 TriaRef::GetParameterValue(&basalfriction,&basalfriction_list[0],gauss_coord); // TO BE DELETED
     5687                TriaRef::GetParameterValue(&basalfriction,&basalfriction_list[0],gauss); // TO BE DELETED
    57805688               
    57815689                /*Calculate scalar parameter*/
    5782                 scalar=gauss_weight*Jdet*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
     5690                scalar=gauss->weight*Jdet*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
    57835691                if(dt){
    57845692                        scalar=dt*scalar;
    … …  
    57945702        VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
    57955703
    5796         /*Free ressources:*/
    5797         xfree((void**)&first_gauss_area_coord);
    5798         xfree((void**)&second_gauss_area_coord);
    5799         xfree((void**)&third_gauss_area_coord);
    5800         xfree((void**)&gauss_weights);
     5704        /*Clean up and return*/
     5705        delete gauss;
     5706        delete friction;
    58015707        xfree((void**)&doflist);
    5802         delete friction;
    58035708}
    58045709/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.