Changeset 6212
- Timestamp:
- 10/08/10 16:34:31 (14 years ago)
- Location:
- issm/trunk/src/c/objects/Elements
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified issm/trunk/src/c/objects/Elements/Penta.cpp ¶
r6200 r6212 2858 2858 2859 2859 /*Intermediaries */ 2860 boolartdiff;2860 int artdiff; 2861 2861 int i,j,ig,found=0; 2862 2862 double Jdet,u,v,w,epsvel; 2863 2863 double gravity,rho_ice,rho_water; 2864 2864 double heatcapacity,thermalconductivity,dt; 2865 double scalar_artdiff; 2866 double tau_parameter,diameter,normu; 2865 2867 double xyz_list[NUMVERTICES][3]; 2866 2868 double B[3][numdof]; … … 2871 2873 double Bprime_advec[3][numdof]; 2872 2874 double L[numdof]; 2875 double dh1dh6[3][6]; 2873 2876 double D_scalar; 2874 2877 double D[3][3]; … … 2903 2906 Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); 2904 2907 Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); 2908 if (artdiff==2) diameter=MinEdgeLength(xyz_list); 2905 2909 2906 2910 /* Start looping on the number of gaussian points: */ … … 2961 2965 /*Artifficial diffusivity*/ 2962 2966 2963 if(artdiff ){2967 if(artdiff==1){ 2964 2968 /*Build K: */ 2965 2969 D_scalar=gauss->weight*Jdet/(pow(u,2)+pow(v,2)+epsvel); … … 2972 2976 MatrixMultiply(&B_artdiff[0][0],2,numdof,1,&K[0][0],2,2,0,&tBD_artdiff[0][0],0); 2973 2977 MatrixMultiply(&tBD_artdiff[0][0],numdof,2,0,&B_artdiff[0][0],2,numdof,0,&Ke_gaussian_artdiff[0][0],0); 2978 } 2979 else if(artdiff==2){ 2980 2981 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss); 2982 2983 normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5); 2984 if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){ 2985 tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity)); 2986 } 2987 else tau_parameter=diameter/(2*normu); 2988 scalar_artdiff=tau_parameter*gauss->weight*Jdet; 2989 2990 for(i=0;i<numdof;i++){ 2991 for(j=0;j<numdof;j++){ 2992 Ke_gaussian_artdiff[i][j]=scalar_artdiff*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i])*(u*dh1dh6[0][j]+v*dh1dh6[1][j]+w*dh1dh6[2][j]); 2993 } 2994 } 2974 2995 } 2975 2996 else{ … … 3746 3767 /*Intermediaries*/ 3747 3768 int i,j,ig,found=0; 3748 int friction_type ;3769 int friction_type,artdiff; 3749 3770 double Jdet,phi,dt; 3750 3771 double rho_ice,heatcapacity; 3772 double thermalconductivity; 3751 3773 double viscosity,temperature; 3774 double tau_parameter,diameter,normu; 3775 double u,v,w,scalar_artdiff; 3752 3776 double scalar_def,scalar_transient; 3753 3777 double temperature_list[NUMVERTICES]; 3754 3778 double xyz_list[NUMVERTICES][3]; 3755 3779 double L[numdof]; 3780 double dh1dh6[3][6]; 3756 3781 double epsilon[6]; 3757 3782 GaussPenta *gauss=NULL; … … 3763 3788 /*Retrieve all inputs and parameters*/ 3764 3789 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); 3765 this->parameters->FindParam(&dt,DtEnum);3766 3790 rho_ice=matpar->GetRhoIce(); 3767 3791 heatcapacity=matpar->GetHeatCapacity(); 3792 thermalconductivity=matpar->GetThermalConductivity(); 3793 this->parameters->FindParam(&dt,DtEnum); 3794 this->parameters->FindParam(&artdiff,ArtDiffEnum); 3768 3795 Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); 3769 3796 Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); … … 3771 3798 Input* temperature_input=NULL; 3772 3799 if (dt) temperature_input=inputs->GetInput(TemperatureEnum); ISSMASSERT(inputs); 3800 if (artdiff==2) diameter=MinEdgeLength(xyz_list); 3773 3801 3774 3802 /* Start looping on the number of gaussian points: */ … … 3795 3823 scalar_transient=temperature*Jdet*gauss->weight; 3796 3824 for(i=0;i<NUMVERTICES;i++) pe->values[i]+=scalar_transient*L[i]; 3825 } 3826 3827 if(artdiff==2){ 3828 GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss); 3829 3830 vx_input->GetParameterValue(&u, gauss); 3831 vy_input->GetParameterValue(&v, gauss); 3832 vz_input->GetParameterValue(&w, gauss); 3833 3834 normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5); 3835 if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){ 3836 tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity)); 3837 } 3838 else tau_parameter=diameter/(2*normu); 3839 3840 scalar_artdiff=tau_parameter*gauss->weight*Jdet*phi/(rho_ice*heatcapacity); 3841 for(i=0;i<NUMVERTICES;i++) pe->values[i]+=scalar_artdiff*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]); 3797 3842 } 3798 3843 } … … 5441 5486 } 5442 5487 /*}}}*/ 5488 /*FUNCTION Penta::MinEdgeLength{{{1*/ 5489 double Penta::MinEdgeLength(double xyz_list[6][3]){ 5490 /*Return the minimum lenght of the nine egdes of the penta*/ 5491 5492 int i,node0,node1; 5493 int edges[9][2]={{0,1},{0,2},{1,2},{3,4},{3,5},{4,5},{0,3},{1,4},{2,5}}; //list of the nine edges 5494 double minlength=-1; 5495 double length; 5496 5497 for(i=0;i<9;i++){ 5498 /*Find the two nodes for this edge*/ 5499 node0=edges[i][0]; 5500 node1=edges[i][1]; 5501 5502 /*Compute the length of this edge and compare it to the minimal length*/ 5503 length=pow(pow(xyz_list[node0][0]-xyz_list[node1][0],2.0)+pow(xyz_list[node0][1]-xyz_list[node1][1],2.0)+pow(xyz_list[node0][2]-xyz_list[node1][2],2.0),0.5); 5504 if(length<minlength || minlength<0) minlength=length; 5505 } 5506 5507 return minlength; 5508 } 5509 /*}}}*/ 5443 5510 /*FUNCTION Penta::ReduceMatrixStokes {{{1*/ 5444 5511 void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){ … … 5604 5671 *(surface_normal+1)=normal[1]/normal_norm; 5605 5672 *(surface_normal+2)=normal[2]/normal_norm; 5606 5607 } 5608 /*}}}*/ 5673 } 5674 /*}}}*/ -
TabularUnified issm/trunk/src/c/objects/Elements/Penta.h ¶
r6200 r6212 218 218 bool IsOnShelf(void); 219 219 bool IsOnWater(void); 220 void ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp); 220 double MinEdgeLength(double xyz_list[6][3]); 221 void ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp); 221 222 void ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp); 222 223 void SetClone(int* minranks); -
TabularUnified issm/trunk/src/c/objects/Elements/Tria.cpp ¶
r6200 r6212 2415 2415 2416 2416 /*Intermediaries */ 2417 boolartdiff;2417 int artdiff; 2418 2418 int i,j,ig,dim; 2419 2419 double Jdettria ,vx,vy,dvxdx,dvydy; … … 2586 2586 2587 2587 /*Intermediaries */ 2588 boolartdiff;2588 int artdiff; 2589 2589 int i,j,ig,dim; 2590 2590 double nx,ny,norm,Jdettria; … … 3157 3157 3158 3158 /*Intermediaries */ 3159 boolartdiff;3159 int artdiff; 3160 3160 int i,j,ig,dim; 3161 3161 double Jdettria,DL_scalar,dt;
Note:
See TracChangeset
for help on using the changeset viewer.