![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 26501
- Timestamp:
- 10/26/21 01:06:03 (3 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 4 edited
-
analyses/StressbalanceAnalysis.cpp (modified) (15 diffs)
-
classes/Elements/Element.h (modified) (1 diff)
-
classes/Elements/Penta.cpp (modified) (1 diff)
-
classes/Elements/Penta.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
r26468 r26501 788 788 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum,0.); 789 789 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum,0.); 790 /*MLHO*/ 790 791 if(isMLHO){ 791 /*itapopo FIXME applying the same initialization for shear vx and shear vy for now*/792 792 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxShearEnum,0.); 793 793 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyShearEnum,0.); 794 /*3D MLHO also need to have VxBase and VyBase for reconstruting Vx and Vy*/ 795 if (iomodel->domaintype==Domain3DEnum) { 796 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxBaseEnum,0.); 797 iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyBaseEnum,0.); 798 } 794 799 } 795 800 if(iomodel->domaintype==Domain2DhorizontalEnum){ … … 2754 2759 ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHO(Element* element){/*{{{*/ 2755 2760 2756 //_error_("Mono Layer Higher-Order called, not fully tested. If you are sure of using it, comment this line.");2757 2758 2761 /* Check if ice in element */ 2759 2762 if(!element->IsIceInElement()) return NULL; 2763 2764 /*Intermediaries*/ 2765 int domaintype; 2766 Element* basalelement; 2767 2768 /*Get basal element*/ 2769 element->FindParam(&domaintype,DomainTypeEnum); 2770 switch(domaintype){ 2771 case Domain2DhorizontalEnum: 2772 basalelement = element; 2773 break; 2774 case Domain3DEnum: 2775 if(!element->IsOnBase()) return NULL; 2776 basalelement = element->SpawnBasalElement(true); 2777 break; 2778 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 2779 } 2780 2760 2781 /*compute all stiffness matrices for this element*/ 2761 ElementMatrix* Ke1=CreateKMatrixMLHOViscous( element);2762 ElementMatrix* Ke2=CreateKMatrixMLHOFriction( element);2782 ElementMatrix* Ke1=CreateKMatrixMLHOViscous(basalelement); 2783 ElementMatrix* Ke2=CreateKMatrixMLHOFriction(basalelement); 2763 2784 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); 2764 2785 2765 2786 /*clean-up and return*/ 2787 if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; 2766 2788 delete Ke1; 2767 2789 delete Ke2; … … 2849 2871 xDelete<IssmDouble>(basis); 2850 2872 return Ke; 2851 2852 //itapopo OLD - testing above2853 /*Intermediaries*/2854 //int domaintype;2855 Element* basalelement;2856 2857 /*Get basal element*/2858 element->FindParam(&domaintype,DomainTypeEnum);2859 switch(domaintype){2860 case Domain2DhorizontalEnum:2861 basalelement = element;2862 break;2863 case Domain3DEnum: case Domain2DverticalEnum:2864 _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");2865 //basalelement = element->SpawnBasalElement();2866 break;2867 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");2868 }2869 2870 //Element* basalelement = element->SpawnBasalElement();2871 //ElementMatrix* Ke = basalelement->NewElementMatrix(MLHOApproximationEnum);2872 ElementMatrix* KeSSA = CreateKMatrixSSAFriction(basalelement); //only to get K11 and K332873 2874 /*Fetch number of nodes and dof for this finite element*/2875 //int numnodes = basalelement->GetNumberOfNodes();2876 2877 for(int i=0;i<numnodes;i++){2878 for(int j=0;j<numnodes;j++){2879 Ke->values[(4*i+0)*2*2*numnodes+4*j+0] = KeSSA->values[2*i*2*numnodes+2*j]; //K112880 Ke->values[(4*i+2)*2*2*numnodes+4*j+2] = KeSSA->values[(2*i+1)*2*numnodes+2*j+1]; //K332881 }2882 }2883 2884 KeSSA->Echo();2885 Ke->Echo();2886 _error_("mesh ");2887 2888 /*Transform Coordinate System*/2889 //basalelement->TransformStiffnessMatrixCoord(Ke,XYMLHOEnum);2890 2891 /*clean-up and return*/2892 if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};2893 delete KeSSA;2894 return Ke;2895 2873 }/*}}}*/ 2896 2874 ElementMatrix* StressbalanceAnalysis::CreateKMatrixMLHOViscous(Element* element){/*{{{*/ 2875 2876 /* Check if ice in element */ 2877 if(!element->IsIceInElement()) return NULL; 2897 2878 2898 2879 /*Intermediaries*/ … … 2900 2881 IssmDouble *xyz_list = NULL; 2901 2882 int domaintype; 2902 Element* basalelement; 2903 2904 /*Get basal element*/ 2905 element->FindParam(&domaintype,DomainTypeEnum); 2906 switch(domaintype){ 2907 case Domain2DhorizontalEnum: 2908 basalelement = element; 2909 break; 2910 case Domain3DEnum: case Domain2DverticalEnum: 2911 _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 2912 //basalelement = element->GetBasalElement()->SpawnBasalElement(true); 2913 break; 2914 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 2915 } 2916 2917 /*Get element on base*/ 2918 //Element* basalelement = element->GetBasalElement()->SpawnBasalElement(true); 2883 int dim=2; 2919 2884 2920 2885 /*Fetch number of nodes and dof for this finite element*/ 2921 int numnodes = basalelement->GetNumberOfNodes();2886 int numnodes = element->GetNumberOfNodes(); 2922 2887 2923 2888 /*Initialize Element matrix and vectors*/ 2924 ElementMatrix* Ke = basalelement->NewElementMatrix(MLHOApproximationEnum);2889 ElementMatrix* Ke = element->NewElementMatrix(MLHOApproximationEnum); 2925 2890 IssmDouble* dbasis = xNew<IssmDouble>(2*numnodes); // like SSA 2926 2891 IssmDouble* basis = xNew<IssmDouble>(numnodes); // like SSA … … 2930 2895 Input* thickness_input= element->GetInput(ThicknessEnum); _assert_(thickness_input); 2931 2896 Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); 2932 //Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); //vertically integrated vx2933 //Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); //vertically integrated vy2934 2897 Input* vxbase_input = element->GetInput(VxBaseEnum); _assert_(vxbase_input); 2935 2898 Input* vybase_input = element->GetInput(VyBaseEnum); _assert_(vybase_input); … … 2940 2903 /* Start looping on the number of gaussian points: */ 2941 2904 Gauss* gauss = element->NewGauss(2); 2942 Gauss* gauss_base = basalelement->NewGauss();2943 2905 while(gauss->next()){ 2944 gauss->SynchronizeGaussBase(gauss_base); 2945 2946 element->JacobianDeterminant(&Jdet,xyz_list,gauss_base); 2947 basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base); 2948 basalelement->NodalFunctions(basis, gauss_base); 2906 2907 element->JacobianDeterminant(&Jdet,xyz_list,gauss); 2908 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); 2909 element->NodalFunctions(basis, gauss); 2949 2910 2950 2911 thickness_input->GetInputValue(&thickness, gauss); 2951 2912 n_input->GetInputValue(&n,gauss); 2952 int dim=2;//itapopo2953 2913 element->material->ViscosityMLHO(&viscosity[0],dim,xyz_list,gauss,vxbase_input,vybase_input,vxshear_input,vyshear_input,thickness_input,n_input); 2954 2914 … … 3021 2981 /*Clean up and return*/ 3022 2982 delete gauss; 3023 delete gauss_base;3024 if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;};3025 2983 xDelete<IssmDouble>(xyz_list); 3026 2984 xDelete<IssmDouble>(dbasis); … … 3196 3154 IssmDouble* surface = xNew<IssmDouble>(numvertices); 3197 3155 3198 element->FindParam(&dim,DomainDimensionEnum);3199 3156 element->FindParam(&domaintype,DomainTypeEnum); 3200 3157 rho_ice =element->FindParam(MaterialsRhoIceEnum); 3201 3158 g =element->FindParam(ConstantsGEnum); 3202 if(dim==2){ 3203 element->GetInputListOnVertices(thickness,ThicknessEnum); 3204 for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i]; 3205 } 3206 else{ 3207 element->GetVerticesCoordinates(&xyz_list); 3208 element->GetInputListOnVertices(surface,SurfaceEnum); 3209 for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]); 3159 switch(domaintype){ 3160 case Domain2DhorizontalEnum: 3161 element->GetInputListOnVertices(thickness,ThicknessEnum); 3162 for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i]; 3163 dim=2; 3164 break; 3165 case Domain3DEnum: 3166 element->GetVerticesCoordinates(&xyz_list); 3167 element->GetInputListOnVertices(surface,SurfaceEnum); 3168 for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]); 3169 dim=2; 3170 break; 3171 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 3210 3172 } 3211 3173 element->AddInput(PressureEnum,pressure,P1Enum); … … 3217 3179 switch(domaintype){ 3218 3180 case Domain2DhorizontalEnum: 3219 // itapopo FIXME use basalelement or element only here3220 3181 basalelement = element; 3221 3182 break; 3222 3183 case Domain3DEnum: 3223 _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 3224 //if(!element->IsOnBase()){xDelete<IssmDouble>(xyz_list); return;} 3225 //basalelement=element->SpawnBasalElement(); 3226 //break; 3184 if(!element->IsOnBase()){xDelete<IssmDouble>(xyz_list); return;} 3185 basalelement=element->SpawnBasalElement(); 3186 break; 3227 3187 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 3228 3188 } … … 3230 3190 /*Fetch number of nodes and dof for this finite element*/ 3231 3191 int numnodes = basalelement->GetNumberOfNodes(); 3232 int numdof = numnodes* 4; //4DOFs per node3192 int numdof = numnodes*dim*2; //2xdim DOFs per node 3233 3193 3234 3194 /*Fetch dof list and allocate solution vectors*/ 3235 3195 basalelement->GetDofListLocal(&doflist,MLHOApproximationEnum,GsetEnum); 3236 //basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);3237 3196 IssmDouble* values = xNew<IssmDouble>(numdof); 3238 3197 IssmDouble* vbx = xNew<IssmDouble>(numnodes); … … 3247 3206 IssmDouble* vel = xNew<IssmDouble>(numnodes); 3248 3207 IssmDouble* n = xNew<IssmDouble>(numnodes); 3208 IssmDouble* H = xNew<IssmDouble>(numnodes); 3209 IssmDouble* s = xNew<IssmDouble>(numnodes); 3249 3210 3250 3211 /*Use the dof list to index into the solution vector: */ 3251 3212 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]]; 3252 3213 3253 //std::cout<<"MLHO ******** Element ID="<<element->Id()<<"\n";3254 //for(i=0;i<numdof;i++){3255 // std::cout<<values[i]*31536000<<"\n";3256 //}3257 3258 3214 /*Transform solution in Cartesian Space*/ 3259 basalelement->TransformSolutionCoord(&values[0],XYEnum); 3260 //basalelement->FindParam(&domaintype,DomainTypeEnum); 3215 if(dim==2) basalelement->TransformSolutionCoord(&values[0],XYEnum); 3261 3216 3262 3217 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ … … 3292 3247 /*Compute the vertically averaged velocities on each node*/ 3293 3248 basalelement->GetInputListOnNodes(&n[0],MaterialsRheologyNEnum,0.); 3294 for(i=0;i<numnodes;i++){ //numnodes of the 2D mesh in which the MLHO is written 3295 vx[i]=vbx[i]+vshx[i]*(n[i]+1)/(n[i]+2); 3296 vy[i]=vby[i]+vshy[i]*(n[i]+1)/(n[i]+2); 3297 } 3298 3299 /*Get Vz and compute vel (vertically averaged vel)*/ 3300 basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.); 3301 for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]); 3302 3303 /*Add vx and vy as inputs to the tria element (vertically averaged velocities): */ 3304 element->AddBasalInput(VxEnum,vx,element->GetElementType()); 3305 element->AddBasalInput(VyEnum,vy,element->GetElementType()); 3306 element->AddBasalInput(VelEnum,vel,element->GetElementType()); 3249 3250 /* Reconstruct vx, vy and vz solutions for 3D problem 3251 * Add vx and vy as inputs to the tria element (vertically averaged velocities): */ 3252 3253 switch(domaintype){ 3254 case Domain2DhorizontalEnum: 3255 for(i=0;i<numnodes;i++){ //numnodes of the 2D mesh in which the MLHO is written 3256 vx[i]=vbx[i]+vshx[i]*(n[i]+1)/(n[i]+2); 3257 vy[i]=vby[i]+vshy[i]*(n[i]+1)/(n[i]+2); 3258 vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i]); 3259 } 3260 element->AddBasalInput(VxEnum,vx,element->GetElementType()); 3261 element->AddBasalInput(VyEnum,vy,element->GetElementType()); 3262 element->AddBasalInput(VelEnum,vel,element->GetElementType()); 3263 break; 3264 case Domain3DEnum: 3265 basalelement->GetInputListOnNodes(&H[0],ThicknessEnum,0.); 3266 basalelement->GetInputListOnNodes(&s[0],SurfaceEnum,0.); 3267 element->Recover3DMLHOInput(VxEnum, numnodes, vbx, vshx, n, H, s); 3268 element->Recover3DMLHOInput(VyEnum, numnodes, vby, vshy, n, H, s); 3269 break; 3270 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 3271 } 3307 3272 3308 3273 /*Free ressources:*/ … … 3320 3285 xDelete<IssmDouble>(xyz_list); 3321 3286 xDelete<IssmDouble>(n); 3287 xDelete<IssmDouble>(s); 3288 xDelete<IssmDouble>(H); 3322 3289 xDelete<int>(doflist); 3323 3290 if(basalelement->IsSpawnedElement()){basalelement->DeleteMaterials(); delete basalelement;}; … … 3333 3300 switch(domaintype){ 3334 3301 case Domain2DhorizontalEnum: dim = 2; dofpernode = 4; break; 3335 case Domain2DverticalEnum: case Domain3DEnum: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); break; 3302 case Domain3DEnum: dim = 2; dofpernode = 4; break; 3303 case Domain2DverticalEnum: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); break; 3336 3304 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet"); 3337 3305 } … … 3361 3329 vxbase_input->GetInputValue(&vbx,gauss); //base vx 3362 3330 vxshear_input->GetInputValue(&vshx,gauss);//shear vx 3363 values[i* 4+0]=vbx; //base vx3364 values[i* 4+1]=vshx; //shear vx3331 values[i*dofpernode+0]=vbx; //base vx 3332 values[i*dofpernode+1]=vshx; //shear vx 3365 3333 vybase_input->GetInputValue(&vby,gauss); //base vy 3366 3334 vyshear_input->GetInputValue(&vshy,gauss);//shear vy 3367 values[i* 4+2]=vby; //base vy3368 values[i* 4+3]=vshy;//shear vy3335 values[i*dofpernode+2]=vby; //base vy 3336 values[i*dofpernode+3]=vshy;//shear vy 3369 3337 } 3370 3338 -
issm/trunk-jpl/src/c/classes/Elements/Element.h
r26495 r26501 343 343 virtual void PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding)=0; 344 344 virtual int PressureInterpolation()=0; 345 virtual void Recover3DMLHOInput(int targetVel_enum, int numnodes, IssmDouble* vb, IssmDouble* vsh, IssmDouble* n, IssmDouble* H, IssmDouble* s){_error_("not implemented yet");}; 345 346 virtual void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe)=0; 346 347 virtual void ResetFSBasalBoundaryCondition()=0; -
issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
r26487 r26501 179 179 else _error_("not implemented yet"); 180 180 } 181 182 } 183 /*}}}*/ 184 void Penta::Recover3DMLHOInput(int targetVel_enum, int numnodes, IssmDouble* vb, IssmDouble* vsh, IssmDouble* n, IssmDouble* H, IssmDouble* s){/*{{{*/ 185 /* Recover the velocity acording to v=vb+(1-\zeta^{n+1})vsh, where \zeta=(s-z)/H 186 * The variables vb, vsh, n, H and s are all from the 2D horizontal mesh(Tria), with "numnodes" DOFs 187 * To project to penta the DOFs are doubled in size 188 * 189 */ 190 _assert_(this->inputs); 191 if(!IsOnBase()) return; 192 else{ 193 if(targetVel_enum==VxEnum || targetVel_enum==VyEnum){ 194 IssmDouble vel[NUMVERTICES]; 195 IssmDouble* xyz_list; 196 Penta* penta = this; 197 _assert_(NUMVERTICES-2*numnodes==0); 198 199 for(;;){ 200 penta->GetVerticesCoordinates(&xyz_list); 201 for(int i=0;i<numnodes;i++) { 202 vel[i] = vb[i] + vsh[i]*(1.0-pow((s[i]-xyz_list[i*3+2])/H[i], (n[i]+1))); 203 vel[i+NUMVERTICES2D] = vb[i] + vsh[i]*(1-pow((s[i]-xyz_list[(i+NUMVERTICES2D)*3+2])/H[i], (n[i]+1))); 204 } 205 206 /*Add to the bottom side of the element*/ 207 penta->AddInput(targetVel_enum,&vel[0],P1Enum); 208 if (penta->IsOnSurface()) break; 209 penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id); 210 } 211 212 xDelete<IssmDouble>(xyz_list); 213 } 214 else _error_("not implemented yet"); 215 } 181 216 182 217 } -
issm/trunk-jpl/src/c/classes/Elements/Penta.h
r26487 r26501 168 168 void PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding); 169 169 int PressureInterpolation(); 170 void Recover3DMLHOInput(int targetVel_enum, int numnodes, IssmDouble* vb, IssmDouble* vsh, IssmDouble* n, IssmDouble* H, IssmDouble* s); 170 171 void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe); 171 172 void ResetFSBasalBoundaryCondition(void);
Note:
See TracChangeset
for help on using the changeset viewer.
