![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 26155
- Timestamp:
- 03/26/21 13:23:35 (4 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 8 edited
-
analyses/StressbalanceAnalysis.cpp (modified) (20 diffs)
-
analyses/StressbalanceAnalysis.h (modified) (1 diff)
-
classes/gauss/GaussTria.cpp (modified) (1 diff)
-
cores/stressbalance_core.cpp (modified) (4 diffs)
-
shared/Enum/Enum.vim (modified) (10 diffs)
-
shared/Enum/EnumDefinitions.h (modified) (1 diff)
-
shared/Enum/EnumToStringx.cpp (modified) (1 diff)
-
shared/Enum/StringToEnumx.cpp (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
r26148 r26155 188 188 } 189 189 else{ 190 IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0); 191 if(iomodel->domaintype!=Domain2DverticalEnum){ 192 IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,1); 190 if(!isMLHO){ 191 IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0); 192 if(iomodel->domaintype!=Domain2DverticalEnum){ 193 IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,1); 194 } 195 } 196 else{//itapopo testing here 197 IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0); 198 IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,1); 199 if(iomodel->domaintype!=Domain2DverticalEnum){ 200 IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2); 201 IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,3); 202 } 193 203 } 194 204 } … … 773 783 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum,0.); 774 784 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum,0.); 785 if(isMLHO){//itapopo 786 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxBaseEnum,0.); 787 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyBaseEnum,0.); 788 } 775 789 iomodel->FetchDataToInput(inputs,elements,"md.stressbalance.loadingforcex",LoadingforceXEnum); 776 790 iomodel->FetchDataToInput(inputs,elements,"md.stressbalance.loadingforcey",LoadingforceYEnum); … … 1211 1225 GetSolutionFromInputsFS(solution,element); 1212 1226 return; 1213 case SSAApproximationEnum: case HOApproximationEnum: case L1L2ApproximationEnum: case MLHOApproximationEnum: caseSIAApproximationEnum:1227 case SSAApproximationEnum: case HOApproximationEnum: case L1L2ApproximationEnum: case SIAApproximationEnum: 1214 1228 GetSolutionFromInputsHoriz(solution,element); 1229 return; 1230 case MLHOApproximationEnum: 1231 GetSolutionFromInputsMLHO(solution,element); 1215 1232 return; 1216 1233 case SSAHOApproximationEnum: case HOFSApproximationEnum: case SSAFSApproximationEnum: … … 2707 2724 /*MLHO*/ 2708 2725 ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/ 2709 _error_("not implemented yet"); 2726 2710 2727 /* Check if ice in element */ 2711 2728 if(!element->IsIceInElement()) return NULL; 2712 2713 2729 /*compute all stiffness matrices for this element*/ 2714 2730 ElementMatrix* Ke1=CreateKMatrixMLHOViscous(element); … … 2724 2740 2725 2741 if(!element->IsOnBase() || element->IsFloating()) return NULL; 2726 Element* basalelement = element->SpawnBasalElement(); 2742 2743 /*Intermediaries*/ 2744 int domaintype; 2745 Element* basalelement; 2746 2747 /*Get basal element*/ 2748 element->FindParam(&domaintype,DomainTypeEnum); 2749 switch(domaintype){ 2750 case Domain2DhorizontalEnum: 2751 basalelement = element; 2752 break; 2753 case Domain3DEnum: case Domain2DverticalEnum: 2754 _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 2755 //basalelement = element->SpawnBasalElement(); 2756 break; 2757 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 2758 } 2759 2760 //Element* basalelement = element->SpawnBasalElement(); 2727 2761 ElementMatrix* Ke = basalelement->NewElementMatrix(MLHOApproximationEnum); 2728 2762 ElementMatrix* KeSSA = CreateKMatrixSSAFriction(basalelement); //only to get K11 and K33 … … 2733 2767 for(int i=0;i<numnodes;i++){ 2734 2768 for(int j=0;j<numnodes;j++){ 2735 Ke->values[( i+0*numnodes)*2*2*numnodes+j+0*numnodes] = KeSSA->values[2*i*2*numnodes+2*j]; //K112736 Ke->values[( i+2*numnodes)*2*2*numnodes+j+2*numnodes] = KeSSA->values[(2*i+1)*2*numnodes+2*j+1]; //K332769 Ke->values[(4*i+0)*2*2*numnodes+4*j+0] = KeSSA->values[2*i*2*numnodes+2*j]; //K11 2770 Ke->values[(4*i+2)*2*2*numnodes+4*j+2] = KeSSA->values[(2*i+1)*2*numnodes+2*j+1]; //K33 2737 2771 } 2738 2772 } 2773 2774 /*Transform Coordinate System*/ 2775 //basalelement->TransformStiffnessMatrixCoord(Ke,XYMLHOEnum); 2739 2776 2740 2777 /*clean-up and return*/ … … 2748 2785 IssmDouble viscosity,Jdet,thickness,n; 2749 2786 IssmDouble *xyz_list = NULL; 2787 int domaintype; 2788 Element* basalelement; 2789 2790 /*Get basal element*/ 2791 element->FindParam(&domaintype,DomainTypeEnum); 2792 switch(domaintype){ 2793 case Domain2DhorizontalEnum: 2794 basalelement = element; 2795 break; 2796 case Domain3DEnum: case Domain2DverticalEnum: 2797 _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 2798 //basalelement = element->GetBasalElement()->SpawnBasalElement(true); 2799 break; 2800 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 2801 } 2750 2802 2751 2803 /*Get element on base*/ 2752 Element* basalelement = element->GetBasalElement()->SpawnBasalElement(true);2804 //Element* basalelement = element->GetBasalElement()->SpawnBasalElement(true); 2753 2805 2754 2806 /*Fetch number of nodes and dof for this finite element*/ … … 2781 2833 n_input->GetInputValue(&n,gauss); 2782 2834 element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input); 2783 2835 //viscosity=10e10;//itapopo 2836 2784 2837 for(int i=0;i<numnodes;i++){//shape functions on tria element 2785 2838 for(int j=0;j<numnodes;j++){ 2786 Ke->values[( i+0*numnodes)*2*2*numnodes+j+0*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2839 Ke->values[(4*i+0)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity*thickness*( 2787 2840 4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] 2788 2841 );//K11 2789 Ke->values[( i+0*numnodes)*2*2*numnodes+j+1*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2842 Ke->values[(4*i+0)*2*2*numnodes+4*j+1] += gauss->weight*Jdet*viscosity*thickness*( 2790 2843 4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] 2791 2844 )*(n+1)/(n+2);//K12 2792 Ke->values[( i+0*numnodes)*2*2*numnodes+j+2*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2845 Ke->values[(4*i+0)*2*2*numnodes+4*j+2] += gauss->weight*Jdet*viscosity*thickness*( 2793 2846 2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i] 2794 2847 );//K13 2795 Ke->values[ (i+0*numnodes)*2*2*numnodes+j+3*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2848 Ke->values[ (4*i+0)*2*2*numnodes+4*j+3] += gauss->weight*Jdet*viscosity*thickness*( 2796 2849 2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i] 2797 2850 )*(n+1)/(n+2);//K14 2798 2851 2799 Ke->values[( i+1*numnodes)*2*2*numnodes+j+0*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2852 Ke->values[(4*i+1)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity*thickness*( 2800 2853 4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] 2801 2854 )*(n+1)/(n+2);//K21 2802 Ke->values[( i+1*numnodes)*2*2*numnodes+j+1*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2855 Ke->values[(4*i+1)*2*2*numnodes+4*j+1] += gauss->weight*Jdet*viscosity*thickness*( 2803 2856 4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] 2804 2857 )*2*pow(n+1,2)/((2*n+3)*(n+2)) … … 2806 2859 gauss->weight*Jdet*viscosity*basis[j]*basis[i]*pow(n+1,2)/(thickness*(2*n+1)) 2807 2860 ;//K22 2808 Ke->values[( i+1*numnodes)*2*2*numnodes+j+2*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2861 Ke->values[(4*i+1)*2*2*numnodes+4*j+2] += gauss->weight*Jdet*viscosity*thickness*( 2809 2862 2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i] 2810 2863 )*(n+1)/(n+2);//K23 2811 Ke->values[( i+1*numnodes)*2*2*numnodes+j+3*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2864 Ke->values[(4*i+1)*2*2*numnodes+4*j+3] += gauss->weight*Jdet*viscosity*thickness*( 2812 2865 2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i] 2813 2866 )*2*pow(n+1,2)/((2*n+3)*(n+2));//K24 2814 2867 2815 Ke->values[( i+2*numnodes)*2*2*numnodes+j+0*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2868 Ke->values[(4*i+2)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity*thickness*( 2816 2869 2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i] 2817 2870 );//K31 2818 Ke->values[( i+2*numnodes)*2*2*numnodes+j+1*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2871 Ke->values[(4*i+2)*2*2*numnodes+4*j+1] += gauss->weight*Jdet*viscosity*thickness*( 2819 2872 2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i] 2820 2873 )*(n+1)/(n+2);//K32 2821 Ke->values[( i+2*numnodes)*2*2*numnodes+j+2*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2874 Ke->values[(4*i+2)*2*2*numnodes+4*j+2] += gauss->weight*Jdet*viscosity*thickness*( 2822 2875 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i] 2823 2876 );//K33 2824 Ke->values[( i+2*numnodes)*2*2*numnodes+j+3*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2877 Ke->values[(4*i+2)*2*2*numnodes+4*j+3] += gauss->weight*Jdet*viscosity*thickness*( 2825 2878 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i] 2826 2879 )*(n+1)/(n+2);//K34 2827 2880 2828 Ke->values[( i+3*numnodes)*2*2*numnodes+j+0*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2881 Ke->values[(4*i+3)*2*2*numnodes+4*j+0] += gauss->weight*Jdet*viscosity*thickness*( 2829 2882 2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i] 2830 2883 )*(n+1)/(n+2);//K41 2831 Ke->values[( i+3*numnodes)*2*2*numnodes+j+1*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2884 Ke->values[(4*i+3)*2*2*numnodes+4*j+1] += gauss->weight*Jdet*viscosity*thickness*( 2832 2885 2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i] 2833 2886 )*2*pow(n+1,2)/((2*n+3)*(n+2));//K42 2834 Ke->values[( i+3*numnodes)*2*2*numnodes+j+2*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2887 Ke->values[(4*i+3)*2*2*numnodes+4*j+2] += gauss->weight*Jdet*viscosity*thickness*( 2835 2888 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i] 2836 2889 )*(n+1)/(n+2);//K43 2837 Ke->values[( i+3*numnodes)*2*2*numnodes+j+3*numnodes] += gauss->weight*Jdet*viscosity*thickness*(2890 Ke->values[(4*i+3)*2*2*numnodes+4*j+3] += gauss->weight*Jdet*viscosity*thickness*( 2838 2891 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i] 2839 2892 )*2*pow(n+1,2)/((2*n+3)*(n+2)) … … 2844 2897 } 2845 2898 } 2846 2899 2847 2900 /*Transform Coordinate System*/ 2848 basalelement->TransformStiffnessMatrixCoord(Ke,XYEnum);2901 //basalelement->TransformStiffnessMatrixCoord(Ke,XYMLHOEnum); 2849 2902 2850 2903 /*Clean up and return*/ … … 2878 2931 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 2879 2932 } 2880 2933 2881 2934 /*compute all load vectors for this element*/ 2882 2935 ElementVector* pe1=CreatePVectorMLHODrivingStress(basalelement); … … 2920 2973 n_input->GetInputValue(&n,gauss); 2921 2974 2922 for(int i=0;i<numnodes;i++){ 2923 pe->values[i +0*numnodes]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]; //F12924 pe->values[i +1*numnodes]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]*(n+1)/(n+2); //F22925 pe->values[i +2*numnodes]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i]; //F32926 pe->values[i +3*numnodes]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i]*(n+1)/(n+2); //F42975 for(int i=0;i<numnodes;i++){// per node: vx (basal vx), vshx, vy (basal vy), vshy 2976 pe->values[i*4+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]; //F1 2977 pe->values[i*4+1]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i]*(n+1)/(n+2); //F2 2978 pe->values[i*4+2]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i]; //F3 2979 pe->values[i*4+3]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i]*(n+1)/(n+2); //F4 2927 2980 } 2928 2981 } 2929 2982 2930 2983 /*Transform coordinate system*/ 2931 element->TransformLoadVectorCoord(pe,XYEnum);2984 //element->TransformLoadVectorCoord(pe,XYMLHOEnum); 2932 2985 2933 2986 /*Clean up and return*/ … … 2993 3046 2994 3047 for (int i=0;i<numnodes;i++){ 2995 pe->values[i +0*numnodes]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; // F12996 pe->values[i +1*numnodes]+= pressure_sh*Jdet*gauss->weight*normal[0]*basis[i]; // F22997 pe->values[i +2*numnodes]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; // F32998 pe->values[i +3*numnodes]+= pressure_sh*Jdet*gauss->weight*normal[1]*basis[i]; // F43048 pe->values[i*4+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i]; // F1 3049 pe->values[i*4+1]+= pressure_sh*Jdet*gauss->weight*normal[0]*basis[i]; // F2 3050 pe->values[i*4+2]+= pressure*Jdet*gauss->weight*normal[1]*basis[i]; // F3 3051 pe->values[i*4+3]+= pressure_sh*Jdet*gauss->weight*normal[1]*basis[i]; // F4 2999 3052 } 3000 3053 } 3001 3054 3002 3055 /*Transform coordinate system*/ 3003 element->TransformLoadVectorCoord(pe,XYEnum);3056 //element->TransformLoadVectorCoord(pe,XYMLHOEnum); 3004 3057 3005 3058 /*Clean up and return*/ … … 3060 3113 3061 3114 /*Fetch dof list and allocate solution vectors*/ 3062 basalelement->GetDofListLocal(&doflist,MLHOApproximationEnum,GsetEnum); //itapopo check this3115 basalelement->GetDofListLocal(&doflist,MLHOApproximationEnum,GsetEnum); 3063 3116 //basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum); 3064 3117 IssmDouble* values = xNew<IssmDouble>(numdof); … … 3073 3126 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]]; 3074 3127 3128 std::cout<<"******** Element ID="<<element->Id()<<"\n"; 3129 for(i=0;i<numdof;i++){ 3130 std::cout<<values[i]<<"\n"; 3131 } 3132 3075 3133 /*Transform solution in Cartesian Space*/ 3076 3134 basalelement->TransformSolutionCoord(&values[0],XYEnum); … … 3079 3137 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ 3080 3138 for(i=0;i<numnodes;i++){ //numnodes of the 2D mesh in which the MLHO is written 3081 vx[i] =values[i +0*numnodes]; //basal vx3082 vshx[i]=values[i +1*numnodes];3139 vx[i] =values[i*4+0]; //basal vx 3140 vshx[i]=values[i*4+1]; 3083 3141 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector"); 3084 3142 if(xIsInf<IssmDouble>(vx[i])) _error_("Inf found in solution vector"); … … 3086 3144 if(xIsInf<IssmDouble>(vshx[i])) _error_("Inf found in solution vector"); 3087 3145 //if(dim==3){ 3088 vy[i] =values[i +2*numnodes];3089 vshy[i]=values[i +3*numnodes];3146 vy[i] =values[i*4+2]; //basal vy 3147 vshy[i]=values[i*4+3]; 3090 3148 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector"); 3091 3149 if(xIsInf<IssmDouble>(vy[i])) _error_("Inf found in solution vector"); … … 3095 3153 } 3096 3154 3155 /*Add vx and vy as inputs to the tria element (basal velocities): */ 3156 element->AddBasalInput(VxBaseEnum,vx,element->GetElementType()); 3157 element->AddBasalInput(VyBaseEnum,vy,element->GetElementType()); 3158 3097 3159 /*Compute suface velocities vx and vy*/ 3098 3160 if(domaintype==Domain2DhorizontalEnum) { … … 3101 3163 } 3102 3164 3103 /*Get Vz and compute vel */3165 /*Get Vz and compute vel (surface)*/ 3104 3166 basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.); 3105 3167 for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]); 3106 3168 3107 /*Add vx and vy as inputs to the tria element : */3169 /*Add vx and vy as inputs to the tria element surface velocities): */ 3108 3170 element->AddBasalInput(VxEnum,vx,element->GetElementType()); 3109 3171 element->AddBasalInput(VyEnum,vy,element->GetElementType()); … … 3121 3183 xDelete<int>(doflist); 3122 3184 if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; 3185 }/*}}}*/ 3186 void StressbalanceAnalysis::GetSolutionFromInputsMLHO(Vector<IssmDouble>* solution,Element* element){/*{{{*/ 3187 3188 IssmDouble vx,vy,vbx,vby; 3189 int domaintype,dim,approximation,dofpernode; 3190 int* doflist = NULL; 3191 3192 /*Get some parameters*/ 3193 element->FindParam(&domaintype,DomainTypeEnum); 3194 switch(domaintype){ 3195 case Domain2DhorizontalEnum: dim = 2; dofpernode = 4; break; 3196 case Domain2DverticalEnum: case Domain3DEnum: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); break; 3197 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 3198 } 3199 3200 /*Fetch number of nodes and dof for this finite element*/ 3201 int numnodes = element->GetNumberOfNodes(); 3202 int numdof = numnodes*dofpernode; 3203 element->GetInputValue(&approximation,ApproximationEnum); 3204 if(approximation!=MLHOApproximationEnum) _error_("mesh "<<EnumToStringx(approximation)<<" not supported here"); 3205 3206 /*Fetch dof list and allocate solution vector*/ 3207 element->GetDofList(&doflist,approximation,GsetEnum); 3208 IssmDouble* values = xNew<IssmDouble>(numdof); 3209 3210 /*Get inputs*/ 3211 Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); 3212 Input* vxbase_input =element->GetInput(VxBaseEnum); _assert_(vxbase_input); 3213 Input* vy_input =NULL; 3214 Input* vybase_input =NULL; 3215 if(domaintype!=Domain2DverticalEnum){ 3216 vy_input =element->GetInput(VyEnum); _assert_(vy_input); 3217 vybase_input=element->GetInput(VyBaseEnum); _assert_(vybase_input); 3218 } 3219 3220 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ 3221 Gauss* gauss=element->NewGauss(); 3222 for(int i=0;i<numnodes;i++){ 3223 gauss->GaussNode(element->FiniteElement(),i); 3224 3225 /*Recover vx and vy*/ 3226 vx_input->GetInputValue(&vx,gauss); 3227 vxbase_input->GetInputValue(&vbx,gauss); 3228 values[i*4+0]=vbx; //basal vx 3229 values[i*4+1]=vx-vbx; //shear vx 3230 //if(dofpernode==2){ 3231 vy_input->GetInputValue(&vy,gauss); 3232 vybase_input->GetInputValue(&vby,gauss); 3233 values[i*4+2]=vby; //basal vy 3234 values[i*4+3]=vy-vby; //shear vy 3235 //} 3236 } 3237 3238 solution->SetValues(numdof,doflist,values,INS_VAL); 3239 3240 /*Free ressources:*/ 3241 delete gauss; 3242 xDelete<IssmDouble>(values); 3243 xDelete<int>(doflist); 3123 3244 }/*}}}*/ 3124 3245 -
issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
r26047 r26155 65 65 ElementVector* CreatePVectorMLHODrivingStress(Element* element); 66 66 void InputUpdateFromSolutionMLHO(IssmDouble* solution,Element* element); 67 void GetSolutionFromInputsMLHO(Vector<IssmDouble>* solution,Element* element); 67 68 /*HO*/ 68 69 ElementMatrix* CreateJacobianMatrixHO(Element* element); -
issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp
r25514 r26155 636 636 } /*}}}*/ 637 637 void GaussTria::SynchronizeGaussBase(Gauss* gauss){/*{{{*/ 638 639 _error_("not supported"); 640 } 641 /*}}}*/ 638 //itapopo check this 639 _assert_(gauss->Enum()==GaussTriaEnum); 640 GaussTria* gauss_tria = xDynamicCast<GaussTria*>(gauss); 641 642 gauss_tria->coord1=this->coord1; 643 gauss_tria->coord2=this->coord2; 644 gauss_tria->coord3=this->coord3; 645 gauss_tria->weight=UNDEF; 646 //_error_("not supported"); 647 } 648 /*}}}*/ -
issm/trunk-jpl/src/c/cores/stressbalance_core.cpp
r25945 r26155 19 19 bool dakota_analysis,control_analysis; 20 20 int domaintype; 21 bool isSIA,isSSA,isL1L2,is HO,isFS,isNitsche;21 bool isSIA,isSSA,isL1L2,isMLHO,isHO,isFS,isNitsche; 22 22 bool save_results; 23 23 int solution_type; … … 31 31 femmodel->parameters->FindParam(&isSSA,FlowequationIsSSAEnum); 32 32 femmodel->parameters->FindParam(&isL1L2,FlowequationIsL1L2Enum); 33 femmodel->parameters->FindParam(&isMLHO,FlowequationIsMLHOEnum); 33 34 femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum); 34 35 femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); … … 74 75 75 76 /*Compute stressbalance for SSA L1L2 HO and FS*/ 76 if(isSSA || isL1L2 || is HO || isFS){77 if(isSSA || isL1L2 || isMLHO || isHO || isFS){ 77 78 analysis = new StressbalanceAnalysis(); 78 79 analysis->Core(femmodel); … … 81 82 82 83 /*Compute vertical velocities*/ 83 if (domaintype==Domain3DEnum && (isSIA || isSSA || isL1L2 || is HO)){84 if (domaintype==Domain3DEnum && (isSIA || isSSA || isL1L2 || isMLHO || isHO)){ 84 85 85 86 /*We need basal melt rates for vertical velocity*/ -
issm/trunk-jpl/src/c/shared/Enum/Enum.vim
r26126 r26155 936 936 syn keyword cConstant VelEnum 937 937 syn keyword cConstant VxAverageEnum 938 syn keyword cConstant VxBaseEnum 938 939 syn keyword cConstant VxEnum 939 940 syn keyword cConstant VxMeshEnum 940 941 syn keyword cConstant VxObsEnum 941 942 syn keyword cConstant VyAverageEnum 943 syn keyword cConstant VyBaseEnum 942 944 syn keyword cConstant VyEnum 943 945 syn keyword cConstant VyMeshEnum … … 1450 1452 syn keyword cType Cfsurfacesquare 1451 1453 syn keyword cType Channel 1454 syn keyword cType classes 1452 1455 syn keyword cType Constraint 1453 1456 syn keyword cType Constraints … … 1456 1459 syn keyword cType ControlInput 1457 1460 syn keyword cType Covertree 1461 syn keyword cType DatasetInput 1458 1462 syn keyword cType DataSetParam 1459 syn keyword cType DatasetInput1460 1463 syn keyword cType Definition 1461 1464 syn keyword cType DependentObject … … 1470 1473 syn keyword cType ElementInput 1471 1474 syn keyword cType ElementMatrix 1475 syn keyword cType Elements 1472 1476 syn keyword cType ElementVector 1473 syn keyword cType Elements1474 1477 syn keyword cType ExponentialVariogram 1475 1478 syn keyword cType ExternalResult … … 1478 1481 syn keyword cType Friction 1479 1482 syn keyword cType Gauss 1483 syn keyword cType GaussianVariogram 1484 syn keyword cType gaussobjects 1480 1485 syn keyword cType GaussPenta 1481 1486 syn keyword cType GaussSeg 1482 1487 syn keyword cType GaussTetra 1483 1488 syn keyword cType GaussTria 1484 syn keyword cType GaussianVariogram1485 1489 syn keyword cType GenericExternalResult 1486 1490 syn keyword cType GenericOption … … 1497 1501 syn keyword cType IssmDirectApplicInterface 1498 1502 syn keyword cType IssmParallelDirectApplicInterface 1503 syn keyword cType krigingobjects 1499 1504 syn keyword cType Load 1500 1505 syn keyword cType Loads … … 1507 1512 syn keyword cType Matice 1508 1513 syn keyword cType Matlitho 1514 syn keyword cType matrixobjects 1509 1515 syn keyword cType MatrixParam 1510 1516 syn keyword cType Misfit … … 1519 1525 syn keyword cType Observations 1520 1526 syn keyword cType Option 1527 syn keyword cType Options 1521 1528 syn keyword cType OptionUtilities 1522 syn keyword cType Options1523 1529 syn keyword cType Param 1524 1530 syn keyword cType Parameters … … 1534 1540 syn keyword cType Regionaloutput 1535 1541 syn keyword cType Results 1542 syn keyword cType Riftfront 1536 1543 syn keyword cType RiftStruct 1537 syn keyword cType Riftfront1538 1544 syn keyword cType SealevelMasks 1539 1545 syn keyword cType Seg 1540 1546 syn keyword cType SegInput 1547 syn keyword cType Segment 1541 1548 syn keyword cType SegRef 1542 syn keyword cType Segment1543 1549 syn keyword cType SpcDynamic 1544 1550 syn keyword cType SpcStatic … … 1559 1565 syn keyword cType Vertex 1560 1566 syn keyword cType Vertices 1561 syn keyword cType classes1562 syn keyword cType gaussobjects1563 syn keyword cType krigingobjects1564 syn keyword cType matrixobjects1565 1567 syn keyword cType AdjointBalancethickness2Analysis 1566 1568 syn keyword cType AdjointBalancethicknessAnalysis -
issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
r26126 r26155 933 933 VelEnum, 934 934 VxAverageEnum, 935 VxBaseEnum, 935 936 VxEnum, 936 937 VxMeshEnum, 937 938 VxObsEnum, 938 939 VyAverageEnum, 940 VyBaseEnum, 939 941 VyEnum, 940 942 VyMeshEnum, -
issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
r26126 r26155 938 938 case VelEnum : return "Vel"; 939 939 case VxAverageEnum : return "VxAverage"; 940 case VxBaseEnum : return "VxBase"; 940 941 case VxEnum : return "Vx"; 941 942 case VxMeshEnum : return "VxMesh"; 942 943 case VxObsEnum : return "VxObs"; 943 944 case VyAverageEnum : return "VyAverage"; 945 case VyBaseEnum : return "VyBase"; 944 946 case VyEnum : return "Vy"; 945 947 case VyMeshEnum : return "VyMesh"; -
issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
r26126 r26155 959 959 else if (strcmp(name,"Vel")==0) return VelEnum; 960 960 else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; 961 else if (strcmp(name,"VxBase")==0) return VxBaseEnum; 961 962 else if (strcmp(name,"Vx")==0) return VxEnum; 962 963 else if (strcmp(name,"VxMesh")==0) return VxMeshEnum; 963 964 else if (strcmp(name,"VxObs")==0) return VxObsEnum; 964 965 else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; 966 else if (strcmp(name,"VyBase")==0) return VyBaseEnum; 965 967 else if (strcmp(name,"Vy")==0) return VyEnum; 966 968 else if (strcmp(name,"VyMesh")==0) return VyMeshEnum; … … 996 998 else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum; 997 999 else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum; 998 else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;999 else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;1000 1000 else stage=9; 1001 1001 } 1002 1002 if(stage==9){ 1003 if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum; 1003 if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum; 1004 else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum; 1005 else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum; 1004 1006 else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum; 1005 1007 else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum; … … 1119 1121 else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum; 1120 1122 else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum; 1121 else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;1122 else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;1123 1123 else stage=10; 1124 1124 } 1125 1125 if(stage==10){ 1126 if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum; 1126 if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum; 1127 else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum; 1128 else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum; 1127 1129 else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum; 1128 1130 else if (strcmp(name,"Channel")==0) return ChannelEnum; … … 1242 1244 else if (strcmp(name,"IntParam")==0) return IntParamEnum; 1243 1245 else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum; 1244 else if (strcmp(name,"Inputs")==0) return InputsEnum;1245 else if (strcmp(name,"Internal")==0) return InternalEnum;1246 1246 else stage=11; 1247 1247 } 1248 1248 if(stage==11){ 1249 if (strcmp(name,"Intersect")==0) return IntersectEnum; 1249 if (strcmp(name,"Inputs")==0) return InputsEnum; 1250 else if (strcmp(name,"Internal")==0) return InternalEnum; 1251 else if (strcmp(name,"Intersect")==0) return IntersectEnum; 1250 1252 else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum; 1251 1253 else if (strcmp(name,"J")==0) return JEnum; … … 1365 1367 else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum; 1366 1368 else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; 1367 else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;1368 else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;1369 1369 else stage=12; 1370 1370 } 1371 1371 if(stage==12){ 1372 if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; 1372 if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; 1373 else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; 1374 else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; 1373 1375 else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; 1374 1376 else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum;
Note:
See TracChangeset
for help on using the changeset viewer.
