![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 22573
- Timestamp:
- 03/20/18 10:38:42 (7 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 2 edited
-
analyses/StressbalanceAnalysis.cpp (modified) (22 diffs)
-
classes/FemModel.cpp (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
r22534 r22573 873 873 parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.FSreconditioning",StressbalanceFSreconditioningEnum)); 874 874 parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.shelf_dampening",StressbalanceShelfDampeningEnum)); 875 parameters->AddObject(iomodel->CopyConstantObject("md.stressbalance.viscosity_overshoot",StressbalanceViscosityOvershootEnum));876 875 parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum)); 877 876 … … 1443 1442 /*Intermediaries*/ 1444 1443 int dim,domaintype; 1445 IssmDouble viscosity,newviscosity,oldviscosity; 1446 IssmDouble viscosity_overshoot,thickness,Jdet; 1444 IssmDouble viscosity,thickness,Jdet; 1447 1445 IssmDouble *xyz_list = NULL; 1448 1446 … … 1468 1466 Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); 1469 1467 Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); 1470 Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input);1471 1468 Input* vy_input = NULL; 1472 Input* vyold_input = NULL;1473 1469 if(dim==2){ 1474 1470 vy_input = element->GetInput(VyEnum); _assert_(vy_input); 1475 vyold_input = element->GetInput(VyPicardEnum); _assert_(vyold_input); 1476 } 1477 element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum); 1471 } 1478 1472 1479 1473 /* Start looping on the number of gaussian points: */ … … 1487 1481 thickness_input->GetInputValue(&thickness, gauss); 1488 1482 element->material->ViscositySSA(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); 1489 element->material->ViscositySSA(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);1490 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);1491 1483 1492 1484 if(dim==2){ 1493 1485 for(int i=0;i<numnodes;i++){ 1494 1486 for(int j=0;j<numnodes;j++){ 1495 Ke->values[2*i*2*numnodes+2*j] += gauss->weight*Jdet* newviscosity*thickness*(1487 Ke->values[2*i*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*thickness*( 1496 1488 4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] 1497 1489 ); 1498 Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet* newviscosity*thickness*(1490 Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*thickness*( 1499 1491 2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i] 1500 1492 ); 1501 Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->weight*Jdet* newviscosity*thickness*(1493 Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*thickness*( 1502 1494 2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i] 1503 1495 ); 1504 Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->weight*Jdet* newviscosity*thickness*(1496 Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*thickness*( 1505 1497 dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] 1506 1498 ); … … 1511 1503 for(int i=0;i<numnodes;i++){ 1512 1504 for(int j=0;j<numnodes;j++){ 1513 Ke->values[i*numnodes+j] += gauss->weight*Jdet* newviscosity*thickness*(1505 Ke->values[i*numnodes+j] += gauss->weight*Jdet*viscosity*thickness*( 1514 1506 4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] 1515 1507 ); … … 1908 1900 } 1909 1901 1910 /*Now, we have to move the previous Vx and Vy inputs to old1911 * status, otherwise, we'll wipe them off: */1912 element->InputChangeName(VxEnum,VxPicardEnum);1913 if(dim==2)element->InputChangeName(VyEnum,VyPicardEnum);1914 1915 1902 /*Add vx and vy as inputs to the tria element: */ 1916 1903 element->AddBasalInput(VxEnum,vx,element->GetElementType()); … … 2236 2223 for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]); 2237 2224 2238 /*Now, we have to move the previous Vx and Vy inputs to old2239 * status, otherwise, we'll wipe them off: */2240 element->InputChangeName(VxEnum,VxPicardEnum);2241 element->InputChangeName(VyEnum,VyPicardEnum);2242 2243 2225 /*Add vx and vy as inputs to the tria element: */ 2244 2226 element->AddBasalInput(VxEnum,vx,element->GetElementType()); … … 2428 2410 /*Intermediaries*/ 2429 2411 int dim,bsize; 2430 IssmDouble viscosity,newviscosity,oldviscosity; 2431 IssmDouble viscosity_overshoot,thickness,Jdet; 2412 IssmDouble viscosity,thickness,Jdet; 2432 2413 IssmDouble *xyz_list = NULL; 2433 2414 … … 2446 2427 element->GetVerticesCoordinates(&xyz_list); 2447 2428 Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); 2448 Input* vxold_input = element->GetInput(VxPicardEnum); _assert_(vxold_input);2449 2429 Input* vy_input = NULL; 2450 Input* vyold_input = NULL;2451 2430 if(dim==3){ 2452 2431 vy_input=element->GetInput(VyEnum); _assert_(vy_input); 2453 vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input); 2454 } 2455 element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum); 2432 } 2456 2433 2457 2434 /* Start looping on the number of gaussian points: */ … … 2463 2440 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); 2464 2441 element->material->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); 2465 element->material->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);2466 2467 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);2468 2442 2469 2443 if(dim==3){ 2470 2444 for(int i=0;i<numnodes;i++){ 2471 2445 for(int j=0;j<numnodes;j++){ 2472 Ke->values[2*i*2*numnodes+2*j] += gauss->weight*Jdet* newviscosity*(2446 Ke->values[2*i*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*( 2473 2447 4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i] 2474 2448 ); 2475 Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet* newviscosity*(2449 Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*( 2476 2450 2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i] 2477 2451 ); 2478 Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->weight*Jdet* newviscosity*(2452 Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*( 2479 2453 2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i] 2480 2454 ); 2481 Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->weight*Jdet* newviscosity*(2455 Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*( 2482 2456 dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i] 2483 2457 ); … … 2488 2462 for(int i=0;i<numnodes;i++){ 2489 2463 for(int j=0;j<numnodes;j++){ 2490 Ke->values[i*numnodes+j] += gauss->weight*Jdet* newviscosity*(2464 Ke->values[i*numnodes+j] += gauss->weight*Jdet*viscosity*( 2491 2465 4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] 2492 2466 ); … … 2905 2879 for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i]); 2906 2880 } 2907 2908 /*Now, we have to move the previous Vx and Vy inputs to old2909 * status, otherwise, we'll wipe them off: */2910 element->InputChangeName(VxEnum,VxPicardEnum);2911 if(dim==3)element->InputChangeName(VyEnum,VyPicardEnum);2912 2881 2913 2882 /*Add vx and vy as inputs to the element: */ … … 5085 5054 else for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i]); 5086 5055 5087 /*Now, we have to move the previous inputs to old5088 * status, otherwise, we'll wipe them off: */5089 element->InputChangeName(VxEnum,VxPicardEnum);5090 element->InputChangeName(VyEnum,VyPicardEnum);5091 if(pnumdof>0) element->InputChangeName(PressureEnum,PressurePicardEnum);5092 if(dim==3) element->InputChangeName(VzEnum,VzPicardEnum);5093 5094 5056 /*Add vx and vy as inputs to the tria element: */ 5095 element->AddInput(VxEnum, vx,element->VelocityInterpolation());5096 element->AddInput(VyEnum, vy,element->VelocityInterpolation());5097 element->AddInput(VelEnum, vel,element->VelocityInterpolation());5057 element->AddInput(VxEnum, vx, element->VelocityInterpolation()); 5058 element->AddInput(VyEnum, vy, element->VelocityInterpolation()); 5059 element->AddInput(VelEnum,vel,element->VelocityInterpolation()); 5098 5060 if(pnumdof>0) element->AddInput(PressureEnum,pressure,element->PressureInterpolation()); 5099 5061 if(dim==3) element->AddInput(VzEnum,vz, element->VelocityInterpolation()); … … 5879 5841 /*Intermediaries */ 5880 5842 int i,j; 5881 IssmDouble Jdet,viscosity ,oldviscosity,newviscosity,viscosity_overshoot; //viscosity5843 IssmDouble Jdet,viscosity; 5882 5844 IssmDouble *xyz_list = NULL; 5883 5845 IssmDouble* B = xNew<IssmDouble>(3*numdofp); … … 5909 5871 /* Get node coordinates and dof list: */ 5910 5872 element->GetVerticesCoordinates(&xyz_list); 5911 element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);5912 5873 Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); 5913 5874 Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); 5914 Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input);5915 Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input);5916 5875 5917 5876 /* Start looping on the number of gaussian points: */ … … 5927 5886 this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria); 5928 5887 element->material->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input); 5929 element->material->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input); 5930 5931 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); 5932 D_scalar=2*newviscosity*gauss->weight*Jdet; 5888 5889 D_scalar=2*viscosity*gauss->weight*Jdet; 5933 5890 for (i=0;i<3;i++) D[i][i]=D_scalar; 5934 5891 … … 6051 6008 int i,j,approximation; 6052 6009 int dim=3; 6053 IssmDouble Jdet,viscosity ,oldviscosity,newviscosity,viscosity_overshoot;6010 IssmDouble Jdet,viscosity; 6054 6011 IssmDouble epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ 6055 6012 IssmDouble epsilons[6]; //6 for FS … … 6071 6028 /*Retrieve all inputs and parameters*/ 6072 6029 element->GetVerticesCoordinates(&xyz_list); 6073 element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);6074 6030 Input* vx_input =element->GetInput(VxEnum); _assert_(vx_input); 6075 6031 Input* vy_input =element->GetInput(VyEnum); _assert_(vy_input); 6076 Input* vxold_input=element->GetInput(VxPicardEnum); _assert_(vxold_input);6077 Input* vyold_input=element->GetInput(VyPicardEnum); _assert_(vyold_input);6078 6032 Input* vz_input =element->GetInput(VzEnum); _assert_(vz_input); 6079 6033 … … 6092 6046 if(approximation==SSAHOApproximationEnum){ 6093 6047 element->material->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input); 6094 element->material->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);6095 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);6096 6048 } 6097 6049 else if (approximation==SSAFSApproximationEnum){ 6098 element->material->ViscosityFS(& newviscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);6050 element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input); 6099 6051 } 6100 6052 else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet"); 6101 6053 6102 D_scalar=2* newviscosity*gauss->weight*Jdet;6054 D_scalar=2*viscosity*gauss->weight*Jdet; 6103 6055 for (i=0;i<3;i++) D[i][i]=D_scalar; 6104 6056 … … 7018 6970 } 7019 6971 7020 /*Now, we have to move the previous Vx and Vy inputs to old7021 * status, otherwise, we'll wipe them off: */7022 element->InputChangeName(VxEnum,VxPicardEnum);7023 element->InputChangeName(VyEnum,VyPicardEnum);7024 element->InputChangeName(VzEnum,VzPicardEnum);7025 element->InputChangeName(PressureEnum,PressurePicardEnum);7026 7027 6972 /*Add vx and vy as inputs to element: */ 7028 6973 element->AddInput(VxEnum,vx,P1Enum); … … 7128 7073 } 7129 7074 7130 /*Now, we have to move the previous Vx and Vy inputs to old7131 * status, otherwise, we'll wipe them off: */7132 element->InputChangeName(VxEnum,VxPicardEnum);7133 element->InputChangeName(VyEnum,VyPicardEnum);7134 element->InputChangeName(VzEnum,VzPicardEnum);7135 element->InputChangeName(PressureEnum,PressurePicardEnum);7136 7137 7075 /*Add vx and vy as inputs to element: */ 7138 7076 element->AddInput(VxEnum,vx,P1Enum); … … 7228 7166 for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]); 7229 7167 7230 /*Now, we have to move the previous Vx and Vy inputs to old7231 * status, otherwise, we'll wipe them off: */7232 element->InputChangeName(VxEnum,VxPicardEnum);7233 element->InputChangeName(VyEnum,VyPicardEnum);7234 element->InputChangeName(PressureEnum,PressurePicardEnum);7235 7236 7168 /*Add vx and vy as inputs to element: */ 7237 7169 element->AddInput(VxEnum,vx,P1Enum); -
issm/trunk-jpl/src/c/classes/FemModel.cpp
r22553 r22573 545 545 int found=-1; 546 546 for(int i=0;i<nummodels;i++){ 547 548 547 if (analysis_type_list[i]==configuration_type){ 549 548 found=i; … … 566 565 if(this->parameters->Exist(ToolkitsOptionsStringsEnum)){ 567 566 ToolkitsOptionsFromAnalysis(this->parameters,analysis_type); 568 if(VerboseSolver()) _printf0_(" toolkits Options set for analysis type: " << EnumToStringx(analysis_type) << "\n"); 569 } 570 571 } 572 /*}}}*/ 567 if(VerboseSolver()) _printf0_(" toolkits Options set for analysis: " << EnumToStringx(analysis_type) << "\n"); 568 } 569 570 }/*}}}*/ 573 571 void FemModel::SetCurrentConfiguration(int configuration_type){/*{{{*/ 574 572 this->SetCurrentConfiguration(configuration_type,configuration_type);
Note:
See TracChangeset
for help on using the changeset viewer.
