Changeset 5741
- Timestamp:
- 09/10/10 09:48:23 (15 years ago)
- Location:
- issm/trunk/src/c/objects
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Tria.cpp
r5738 r5741 2605 2605 2606 2606 } 2607 } // for (ig=0; ig<num_gauss; ig++)2607 } 2608 2608 2609 2609 /*Add Ke_gg to global matrix Kgg: */ … … 2692 2692 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg2[i][j]; 2693 2693 2694 } // for (ig=0; ig<num_gauss; ig++)2694 } 2695 2695 2696 2696 /*Add Ke_gg to global matrix Kgg: */ … … 2859 2859 } 2860 2860 2861 } // for (ig=0; ig<num_gauss; ig++)2861 } 2862 2862 2863 2863 /*Add Ke_gg to global matrix Kgg: */ … … 2982 2982 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j]; 2983 2983 2984 } // for (ig=0; ig<num_gauss; ig++)2984 } 2985 2985 2986 2986 /*Add Ke_gg to global matrix Kgg: */ … … 3109 3109 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j]; 3110 3110 3111 } // for (ig=0; ig<num_gauss; ig++)3111 } 3112 3112 3113 3113 /*Add Ke_gg to global matrix Kgg: */ … … 3231 3231 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j]; 3232 3232 3233 } // for (ig=0; ig<num_gauss; ig++)3233 } 3234 3234 3235 3235 /*Add Ke_gg to global matrix Kgg: */ … … 3351 3351 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j]; 3352 3352 3353 } // for (ig=0; ig<num_gauss; ig++)3353 } 3354 3354 3355 3355 /*Add Ke_gg to global matrix Kgg: */ … … 3488 3488 3489 3489 3490 } //for (ig=0; ig<num_gauss; ig++)3490 } 3491 3491 3492 3492 /*Add Ke_gg to global matrix Kgg: */ … … 3725 3725 } 3726 3726 3727 } // for (ig=0; ig<num_gauss; ig++)3727 } 3728 3728 3729 3729 /*Add Ke_gg to global matrix Kgg: */ … … 3835 3835 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg2[i][j]; 3836 3836 3837 } // for (ig=0; ig<num_gauss; ig++)3837 } 3838 3838 3839 3839 /*Add Ke_gg to global matrix Kgg: */ … … 4034 4034 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g-dhdt_g)*L[i]; 4035 4035 4036 } // for (ig=0; ig<num_gauss; ig++)4036 } 4037 4037 4038 4038 /*Add pe_g to global matrix Kgg: */ … … 4104 4104 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g-dhdt_g)*L[i]; 4105 4105 4106 } // for (ig=0; ig<num_gauss; ig++)4106 } 4107 4107 4108 4108 /*Add pe_g to global matrix Kgg: */ … … 4170 4170 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g)*L[i]; 4171 4171 4172 } // for (ig=0; ig<num_gauss; ig++)4172 } 4173 4173 4174 4174 /*Add pe_g to global matrix Kgg: */ … … 4275 4275 void Tria::CreatePVectorDiagnosticMacAyeal( Vec pg){ 4276 4276 4277 int i,j ;4277 int i,j,ig; 4278 4278 4279 4279 /* node data: */ 4280 const int numdof= 2*NUMVERTICES;4280 const int numdof=NDOF2*NUMVERTICES; 4281 4281 double xyz_list[NUMVERTICES][3]; 4282 4282 int* doflist=NULL; … … 4284 4284 /* parameters: */ 4285 4285 double plastic_stress; 4286 double slope[ NDOF2];4286 double slope[2]; 4287 4287 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 4299 4290 double Jdet; 4300 4301 /*nodal functions: */4302 4291 double l1l2l3[3]; 4303 4304 /*element vector at the gaussian points: */4305 4292 double pe_g[numdof]={0.0}; 4306 4293 double pe_g_gaussian[numdof]; 4307 4308 /*input parameters for structural analysis (diagnostic): */4309 4294 double thickness; 4310 4311 /*inputs: */4312 4295 bool onwater; 4313 4296 int drag_type; 4314 Input* thickness_input=NULL;4315 Input* surface_input=NULL;4316 Input* drag_input=NULL;4317 4297 4318 4298 /*retrieve inputs :*/ 4319 4299 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); 4320 4300 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); 4324 4304 4325 4305 /*First, if we are on water, return empty vector: */ 4326 if(onwater) return;4306 if(onwater) return; 4327 4307 4328 4308 /* Get node coordinates and dof list: */ … … 4330 4310 GetDofList(&doflist,MacAyealApproximationEnum); 4331 4311 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 4336 4312 /* 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); 4347 4320 4348 4321 /*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the 4349 4322 * element itself: */ 4350 4323 if(drag_type==1){ 4351 drag_input->GetParameterValue(&plastic_stress, &gauss_l1l2l3[0]);4324 drag_input->GetParameterValue(&plastic_stress,gauss); 4352 4325 } 4353 4326 4354 4327 /* Get Jacobian determinant: */ 4355 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss _l1l2l3);4328 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); 4356 4329 4357 4330 /*Get nodal functions: */ 4358 GetNodalFunctions(l1l2l3, gauss _l1l2l3);4331 GetNodalFunctions(l1l2l3, gauss); 4359 4332 4360 4333 /*Compute driving stress: */ … … 4365 4338 for (i=0;i<NUMVERTICES;i++){ 4366 4339 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]; 4368 4341 } 4369 4342 } … … 4372 4345 for (i=0;i<NUMVERTICES;i++){ 4373 4346 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]; 4375 4348 } 4376 4349 } … … 4380 4353 for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i]; 4381 4354 4382 } //for (ig=0; ig<num_gauss; ig++)4355 } 4383 4356 4384 4357 /*Add pe_g to global vector pg: */ … … 4386 4359 4387 4360 /*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; 4392 4362 xfree((void**)&doflist); 4393 4363 } … … 5043 5013 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i]; 5044 5014 5045 } // for (ig=0; ig<num_gauss; ig++)5015 } 5046 5016 5047 5017 /*Add pe_g to global matrix Kgg: */ … … 5116 5086 for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss->weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i]; 5117 5087 5118 } // for (ig=0; ig<num_gauss; ig++)5088 } 5119 5089 5120 5090 /*Add pe_g to global matrix Kgg: */ … … 5195 5165 for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i]; 5196 5166 5197 } //for (ig=0; ig<num_gauss; ig++)5167 } 5198 5168 5199 5169 /*Add pe_g to global vector pg: */ … … 5766 5736 void Tria::GradjDragStokes(Vec gradient){ 5767 5737 5768 int i ;5738 int i,ig; 5769 5739 5770 5740 /* node data: */ … … 5777 5747 double alpha_complement; 5778 5748 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; 5788 5750 5789 5751 /* parameters: */ … … 5798 5760 double grade_g[NUMVERTICES]={0.0}; 5799 5761 double grade_g_gaussian[NUMVERTICES]; 5800 5801 /* Jacobian: */5802 5762 double Jdet; 5803 5804 /*nodal functions: */5805 5763 double l1l2l3[3]; 5806 5807 /* strain rate: */5808 5764 double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ 5809 5810 /*inputs: */5811 5765 bool shelf; 5812 5766 int drag_type; 5813 5814 /*parameters: */5815 5767 double cm_noisedmp; 5816 5817 5768 int analysis_type; 5818 5769 … … 5823 5774 inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); 5824 5775 inputs->GetParameterValue(&drag_type,DragTypeEnum); 5776 Input* drag_input=inputs->GetInput(DragCoefficientEnum); ISSMASSERT(drag_input); 5825 5777 5826 5778 /*retrieve some parameters: */ … … 5838 5790 friction=new Friction("2d",inputs,matpar,analysis_type); 5839 5791 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 5843 5792 /* 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); 5850 5797 5851 5798 /*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); 5853 5800 else alpha_complement=0; 5854 inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum);5801 inputs->GetParameterValue(&drag,gauss,DragCoefficientEnum); 5855 5802 5856 5803 /*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); 5860 5807 5861 5808 /*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); 5865 5812 5866 5813 /*Get normal vector to the bed */ … … 5872 5819 5873 5820 /* Get Jacobian determinant: */ 5874 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss _l1l2l3);5821 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss); 5875 5822 5876 5823 /* Get nodal functions value at gaussian point:*/ 5877 GetNodalFunctions(l1l2l3, gauss _l1l2l3);5824 GetNodalFunctions(l1l2l3, gauss); 5878 5825 5879 5826 /*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); 5881 5828 5882 5829 /*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); 5884 5831 5885 5832 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ … … 5890 5837 -mu *(2*drag*alpha_complement*(vy - vz*bed_normal[1]*bed_normal[2])) 5891 5838 -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]; 5893 5840 5894 5841 //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]); 5896 5843 } 5897 5844 … … 5903 5850 VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES); 5904 5851 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);5912 5852 delete friction; 5853 delete gauss; 5913 5854 5914 5855 } -
issm/trunk/src/c/objects/Loads/Numericalflux.cpp
r5739 r5741 652 652 int count=0; 653 653 int type; 654 int numberofnodes=2; 655 656 /*output: */ 654 int numberofnodes; 657 655 int* doflist=NULL; 658 656 659 /*recover type: */660 inputs->GetParameterValue(&type,TypeEnum);661 662 657 /*Some checks for debugging*/ 663 658 ISSMASSERT(nodes); 664 659 665 660 /*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 } 669 670 670 671 /*Figure out size of doflist: */ … … 691 692 692 693 /*Build unit outward pointing vector*/ 693 const int numgrids=4;694 694 double vector[2]; 695 695 double norm;
Note:
See TracChangeset
for help on using the changeset viewer.