Changeset 13051
- Timestamp:
- 08/15/12 17:15:45 (13 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/Makefile.am
r13046 r13051 203 203 ./shared/Numerics/isnan.h\ 204 204 ./shared/Numerics/isnan.cpp\ 205 ./shared/Numerics/cubic.cpp\ 205 206 ./shared/Numerics/extrema.cpp\ 206 207 ./shared/Numerics/XZvectorsToCoordinateSystem.cpp\ -
issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp
r13042 r13051 1677 1677 //Need to know the type of approximation for this element 1678 1678 if(iomodel->Data(FlowequationElementEquationEnum)){ 1679 if ( *(iomodel->Data(FlowequationElementEquationEnum)+index)==MacAyealApproximationEnum){1679 if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealApproximationEnum){ 1680 1680 this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealApproximationEnum)); 1681 1681 } 1682 else if ( *(iomodel->Data(FlowequationElementEquationEnum)+index)==PattynApproximationEnum){1682 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==PattynApproximationEnum){ 1683 1683 this->inputs->AddInput(new IntInput(ApproximationEnum,PattynApproximationEnum)); 1684 1684 } 1685 else if ( *(iomodel->Data(FlowequationElementEquationEnum)+index)==MacAyealPattynApproximationEnum){1685 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealPattynApproximationEnum){ 1686 1686 this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealPattynApproximationEnum)); 1687 1687 } 1688 else if ( *(iomodel->Data(FlowequationElementEquationEnum)+index)==HutterApproximationEnum){1688 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==HutterApproximationEnum){ 1689 1689 this->inputs->AddInput(new IntInput(ApproximationEnum,HutterApproximationEnum)); 1690 1690 } 1691 else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==StokesApproximationEnum){ 1691 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==L1L2ApproximationEnum){ 1692 this->inputs->AddInput(new IntInput(ApproximationEnum,L1L2ApproximationEnum)); 1693 } 1694 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==StokesApproximationEnum){ 1692 1695 this->inputs->AddInput(new IntInput(ApproximationEnum,StokesApproximationEnum)); 1693 1696 } 1694 else if ( *(iomodel->Data(FlowequationElementEquationEnum)+index)==MacAyealStokesApproximationEnum){1697 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealStokesApproximationEnum){ 1695 1698 this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealStokesApproximationEnum)); 1696 1699 } 1697 else if ( *(iomodel->Data(FlowequationElementEquationEnum)+index)==PattynStokesApproximationEnum){1700 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==PattynStokesApproximationEnum){ 1698 1701 this->inputs->AddInput(new IntInput(ApproximationEnum,PattynStokesApproximationEnum)); 1699 1702 } 1700 else if ( *(iomodel->Data(FlowequationElementEquationEnum)+index)==NoneApproximationEnum){1703 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==NoneApproximationEnum){ 1701 1704 this->inputs->AddInput(new IntInput(ApproximationEnum,NoneApproximationEnum)); 1702 1705 } 1703 1706 else{ 1704 _error_("Approximation type " << EnumToStringx((int) *(iomodel->Data(FlowequationElementEquationEnum)+index)) << " not supported yet");1707 _error_("Approximation type " << EnumToStringx((int)iomodel->Data(FlowequationElementEquationEnum)[index]) << " not supported yet"); 1705 1708 } 1706 1709 } … … 2801 2804 2802 2805 /*Intermediaries*/ 2803 IssmInt i,j;2804 int penta_type;2805 int penta_node_ids[6];2806 int penta_vertex_ids[6];2807 IssmDouble 2808 IssmDouble 2809 int stabilization;2810 bool dakota_analysis;2811 bool isstokes;2812 IssmDouble 2806 IssmInt i,j; 2807 int penta_type; 2808 int penta_node_ids[6]; 2809 int penta_vertex_ids[6]; 2810 IssmDouble nodeinputs[6]; 2811 IssmDouble yts; 2812 int stabilization; 2813 bool dakota_analysis; 2814 bool isstokes; 2815 IssmDouble beta,heatcapacity,referencetemperature,meltingpoint,latentheat; 2813 2816 2814 2817 /*Fetch parameters: */ … … 4232 4235 const int numdof=NDOF1*NUMVERTICES; 4233 4236 4234 int 4235 int * doflist=NULL;4236 IssmDouble 4237 IssmDouble 4238 GaussPenta *gauss=NULL;4237 int i; 4238 int *doflist = NULL; 4239 IssmDouble values[numdof]; 4240 IssmDouble temp; 4241 GaussPenta *gauss = NULL; 4239 4242 4240 4243 /*Get dof list: */ … … 6200 6203 case MacAyealApproximationEnum: 6201 6204 return CreateKMatrixDiagnosticMacAyeal2d(); 6205 case L1L2ApproximationEnum: 6206 return CreateKMatrixDiagnosticL1L2(); 6202 6207 case PattynApproximationEnum: 6203 6208 return CreateKMatrixDiagnosticPattyn(); … … 6226 6231 6227 6232 /*Intermediaries*/ 6228 int connectivity[2];6229 int i,i0,i1,j0,j1;6230 IssmDouble 6233 int connectivity[2]; 6234 int i,i0,i1,j0,j1; 6235 IssmDouble one0,one1; 6231 6236 6232 6237 /*Initialize Element matrix*/ … … 6453 6458 delete Ke2; 6454 6459 delete Ke3; 6460 return Ke; 6461 } 6462 /*}}}*/ 6463 /*FUNCTION Penta::CreateKMatrixDiagnosticL1L2{{{*/ 6464 ElementMatrix* Penta::CreateKMatrixDiagnosticL1L2(void){ 6465 6466 /*compute all stiffness matrices for this element*/ 6467 ElementMatrix* Ke1=CreateKMatrixDiagnosticL1L2Viscous(); 6468 ElementMatrix* Ke2=CreateKMatrixDiagnosticL1L2Friction(); 6469 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); 6470 6471 /*clean-up and return*/ 6472 delete Ke1; 6473 delete Ke2; 6474 return Ke; 6475 } 6476 /*}}}*/ 6477 /*FUNCTION Penta::CreateKMatrixDiagnosticL1L2Viscous{{{*/ 6478 ElementMatrix* Penta::CreateKMatrixDiagnosticL1L2Viscous(void){ 6479 6480 /*Constants*/ 6481 const int numdof2d=2*NUMVERTICES2D; 6482 6483 /*Intermediaries */ 6484 int i,j; 6485 IssmDouble Jdet,viscosity; 6486 IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ 6487 IssmDouble xyz_list[NUMVERTICES][3]; 6488 IssmDouble B[3][numdof2d]; 6489 IssmDouble Bprime[3][numdof2d]; 6490 IssmDouble Ke_gg_gaussian[numdof2d][numdof2d]; 6491 IssmDouble D[3][3]= {0.0}; // material matrix, simple scalar matrix. 6492 Tria *tria = NULL; 6493 Penta *pentabase = NULL; 6494 GaussPenta *gauss = NULL; 6495 GaussTria *gauss_tria = NULL; 6496 6497 /*Find penta on bed as this is a macayeal elements: */ 6498 pentabase=GetBasalElement(); 6499 tria=pentabase->SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 6500 6501 /*Initialize Element matrix*/ 6502 ElementMatrix* Ke=new ElementMatrix(tria->nodes,NUMVERTICES2D,this->parameters,L1L2ApproximationEnum); 6503 6504 /*Retrieve all inputs and parameters*/ 6505 GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES); 6506 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); 6507 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); 6508 Input* surf_input=inputs->GetInput(SurfaceEnum); _assert_(surf_input); 6509 6510 /* Start looping on the number of gaussian points: */ 6511 gauss=new GaussPenta(5,5); 6512 gauss_tria=new GaussTria(); 6513 for(int ig=gauss->begin();ig<gauss->end();ig++){ 6514 6515 gauss->GaussPoint(ig); 6516 gauss->SynchronizeGaussTria(gauss_tria); 6517 6518 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); 6519 tria->GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss_tria); 6520 tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria); 6521 6522 /*Get viscosity for L1L2 model*/ 6523 GetL1L2Viscosity(&viscosity,&xyz_list[0][0],gauss,vx_input,vy_input,surf_input); 6524 6525 for(i=0;i<3;i++) D[i][i]=2*viscosity*gauss->weight*Jdet; 6526 6527 TripleMultiply( &B[0][0],3,numdof2d,1, 6528 &D[0][0],3,3,0, 6529 &Bprime[0][0],3,numdof2d,0, 6530 &Ke_gg_gaussian[0][0],0); 6531 for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof2d+j]+=Ke_gg_gaussian[i][j]; 6532 } 6533 6534 /*Transform Coordinate System*/ 6535 TransformStiffnessMatrixCoord(Ke,tria->nodes,NUMVERTICES2D,XYEnum); 6536 6537 /*Clean up and return*/ 6538 delete tria->matice; 6539 delete tria; 6540 delete gauss_tria; 6541 delete gauss; 6542 return Ke; 6543 } 6544 /*}}}*/ 6545 /*FUNCTION Penta::CreateKMatrixDiagnosticL1L2Friction{{{*/ 6546 ElementMatrix* Penta::CreateKMatrixDiagnosticL1L2Friction(void){ 6547 6548 /*Initialize Element matrix and return if necessary*/ 6549 if(IsFloating() || !IsOnBed()) return NULL; 6550 6551 /*Build a tria element using the 3 nodes of the base of the penta. Then use 6552 * the tria functionality to build a friction stiffness matrix on these 3 6553 * nodes: */ 6554 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 6555 ElementMatrix* Ke=tria->CreateKMatrixDiagnosticMacAyealFriction(); 6556 delete tria->matice; delete tria; 6557 6558 /*clean-up and return*/ 6455 6559 return Ke; 6456 6560 } … … 7193 7297 case PattynApproximationEnum: 7194 7298 return CreatePVectorDiagnosticPattyn(); 7299 case L1L2ApproximationEnum: 7300 return CreatePVectorDiagnosticL1L2(); 7195 7301 case HutterApproximationEnum: 7196 7302 return NULL; … … 7346 7452 /*FUNCTION Penta::CreatePVectorDiagnosticMacAyeal{{{*/ 7347 7453 ElementVector* Penta::CreatePVectorDiagnosticMacAyeal(void){ 7454 7455 if (!IsOnBed()) return NULL; 7456 7457 /*Call Tria function*/ 7458 Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. 7459 ElementVector* pe=tria->CreatePVectorDiagnosticMacAyeal(); 7460 delete tria->matice; delete tria; 7461 7462 /*Clean up and return*/ 7463 return pe; 7464 } 7465 /*}}}*/ 7466 /*FUNCTION Penta::CreatePVectorDiagnosticL1L2{{{*/ 7467 ElementVector* Penta::CreatePVectorDiagnosticL1L2(void){ 7348 7468 7349 7469 if (!IsOnBed()) return NULL; … … 7892 8012 const int numdof=NDOF2*NUMVERTICES; 7893 8013 7894 int 7895 int 7896 int * doflist=NULL;7897 IssmDouble 7898 IssmDouble 7899 GaussPenta *gauss;8014 int i; 8015 int approximation; 8016 int *doflist = NULL; 8017 IssmDouble vx,vy; 8018 IssmDouble values[numdof]; 8019 GaussPenta *gauss; 7900 8020 7901 8021 /*Get approximation enum and dof list: */ … … 8043 8163 } 8044 8164 /*}}}*/ 8165 /*FUNCTION Penta::GetL1L2Viscosity{{{*/ 8166 void Penta::GetL1L2Viscosity(IssmDouble* pviscosity,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input,Input* surface_input){ 8167 /*Compute the L1L2 viscosity 8168 * 8169 * 1 8170 * mu = - A^-1 (sigma'_e)^(1-n) 8171 * 2 8172 * 8173 * sigma'_e^2 = |sigma'_//|^2 + |sigma'_perp|^2 (see Perego 2012 eq. 17,18) 8174 * 8175 * L1L2 assumptions: 8176 * 8177 * (1) |eps_b|_// = A (|sigma'_//|^2 + |sigma'_perp|^2)^((n-1)/2) |sigma'_//| 8178 * (2) |sigma'_perp|^2 = |rho g (s-z) grad(s)|^2 8179 * 8180 * Assuming that n = 3, we have a polynom of degree 3 to solve (the only unkown is X=|sigma'_//|) 8181 * 8182 * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0 */ 8183 8184 _error_("Not supported yet"); 8185 int i; 8186 IssmDouble a,c,d,z,s,viscosity; 8187 IssmDouble tau_perp,tau_par,eps_b,A; 8188 IssmDouble epsilonvx[5]; /*exx eyy exy exz eyz*/ 8189 IssmDouble epsilonvy[5]; /*exx eyy exy exz eyz*/ 8190 IssmDouble epsilon[5]; /*exx eyy exy exz eyz*/ 8191 IssmDouble z_list[NUMVERTICES]; 8192 IssmDouble slope[3]; 8193 8194 /*Check that both inputs have been found*/ 8195 if (!vx_input || !vy_input || !surface_input) _error_("Input missing"); 8196 8197 /*Get tau_perp*/ 8198 for(i=0;i<NUMVERTICES;i++) z_list[i]=xyz_list[3*i+2]; 8199 surface_input->GetInputValue(&s,gauss); 8200 surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); 8201 PentaRef::GetInputValue(&z,&z_list[0],gauss); 8202 tau_perp = matpar->GetRhoIce() * matpar->GetG() * (s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]); 8203 8204 /* Get eps_b*/ 8205 vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss); 8206 vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss); 8207 for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; 8208 eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]); 8209 if(eps_b==0.){ 8210 *pviscosity=2.5*pow(10.,17.); 8211 return; 8212 } 8213 8214 /*Get A*/ 8215 _assert_(matice->GetN()==3.0); 8216 A=matice->GetA(); 8217 8218 /*Solve for tau_par (http://en.wikipedia.org/wiki/Cubic_function)*/ 8219 a=A; 8220 c=A*tau_perp*tau_perp; 8221 d=-eps_b; 8222 tau_par = -1./(3.*a) * pow(0.5*(27.*a*a*d + sqrt(27.*a*a*d*27.*a*a*d -4.*pow(-3.*a*c,3.))),1./3.); 8223 8224 //printf("=======================================\n"); 8225 //printf("tau_par = %g (%g=0?)\n",tau_par,a*tau_par*tau_par*tau_par+c*tau_par+d); 8226 //printf("a=%g c=%g d=%g\n",a,c,d); 8227 //IssmDouble coeff[4]; 8228 //coeff[0]=d; 8229 //coeff[1]=c; 8230 //coeff[2]=0.; 8231 //coeff[3]=a; 8232 //int numroots; 8233 //IssmDouble roots[3]; 8234 //cubic(coeff,roots,&numroots); 8235 //tau_par=roots[0]; 8236 //for(i=0;i<numroots;i++) printf(" %g ",roots[i]); 8237 //printf(" (%g =0?)\n",a*roots[0]*roots[0]*roots[0]+c*roots[0]+d); 8238 8239 8240 /*Viscosity*/ 8241 viscosity = 1./(2.*A)*pow(tau_par*tau_par + tau_perp*tau_perp ,-2.); 8242 _assert_(!isnan(viscosity)); 8243 8244 /*Assign output pointer*/ 8245 *pviscosity = viscosity; 8246 return; 8247 } 8248 /*}}}*/ 8045 8249 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHoriz {{{*/ 8046 8250 void Penta::InputUpdateFromSolutionDiagnosticHoriz(IssmDouble* solution){ … … 8061 8265 return; 8062 8266 } 8267 } 8268 if (approximation==L1L2ApproximationEnum){ 8269 if (!IsOnBed()) return; 8270 InputUpdateFromSolutionDiagnosticL1L2(solution); 8271 return; 8063 8272 } 8064 8273 else if (approximation==PattynApproximationEnum){ … … 8342 8551 xDelete<int>(doflistm); 8343 8552 xDelete<int>(doflists); 8553 } 8554 /*}}}*/ 8555 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticL1L2 {{{*/ 8556 void Penta::InputUpdateFromSolutionDiagnosticL1L2(IssmDouble* solution){ 8557 8558 const int numdof=NDOF2*NUMVERTICES; 8559 8560 int i; 8561 IssmDouble rho_ice,g; 8562 IssmDouble values[numdof]; 8563 IssmDouble vx[NUMVERTICES]; 8564 IssmDouble vy[NUMVERTICES]; 8565 IssmDouble vz[NUMVERTICES]; 8566 IssmDouble vel[NUMVERTICES]; 8567 IssmDouble pressure[NUMVERTICES]; 8568 IssmDouble surface[NUMVERTICES]; 8569 IssmDouble xyz_list[NUMVERTICES][3]; 8570 int *doflist = NULL; 8571 Penta *penta = NULL; 8572 8573 /*Get dof list: */ 8574 GetDofList(&doflist,L1L2ApproximationEnum,GsetEnum); 8575 8576 /*Use the dof list to index into the solution vector: */ 8577 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]]; 8578 8579 /*Transform solution in Cartesian Space*/ 8580 TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,XYEnum); /*2D: only the first 3 nodes are taken*/ 8581 8582 /*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */ 8583 for(i=0;i<3;i++){ 8584 vx[i] =values[i*NDOF2+0]; 8585 vy[i] =values[i*NDOF2+1]; 8586 vx[i+3]=vx[i]; 8587 vy[i+3]=vy[i]; 8588 8589 /*Check solution*/ 8590 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector"); 8591 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector"); 8592 } 8593 8594 /*Get parameters fro pressure computation*/ 8595 rho_ice=matpar->GetRhoIce(); 8596 g=matpar->GetG(); 8597 8598 /*Start looping over all elements above current element and update all inputs*/ 8599 penta=this; 8600 for(;;){ 8601 8602 /*Get node data: */ 8603 GetVerticesCoordinates(&xyz_list[0][0],penta->nodes,NUMVERTICES); 8604 8605 /*Now Compute vel*/ 8606 GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0 8607 for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5); 8608 8609 /*Now compute pressure*/ 8610 GetInputListOnVertices(&surface[0],SurfaceEnum); 8611 for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]); 8612 8613 /*Now, we have to move the previous Vx and Vy inputs to old 8614 * status, otherwise, we'll wipe them off: */ 8615 penta->inputs->ChangeEnum(VxEnum,VxPicardEnum); 8616 penta->inputs->ChangeEnum(VyEnum,VyPicardEnum); 8617 penta->inputs->ChangeEnum(PressureEnum,PressurePicardEnum); 8618 8619 /*Add vx and vy as inputs to the tria element: */ 8620 penta->inputs->AddInput(new PentaP1Input(VxEnum,vx)); 8621 penta->inputs->AddInput(new PentaP1Input(VyEnum,vy)); 8622 penta->inputs->AddInput(new PentaP1Input(VelEnum,vel)); 8623 penta->inputs->AddInput(new PentaP1Input(PressureEnum,pressure)); 8624 8625 /*Stop if we have reached the surface*/ 8626 if (penta->IsOnSurface()) break; 8627 8628 /* get upper Penta*/ 8629 penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id); 8630 } 8631 8632 /*Free ressources:*/ 8633 xDelete<int>(doflist); 8344 8634 } 8345 8635 /*}}}*/ -
issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h
r13036 r13051 227 227 ElementMatrix* CreateKMatrixDiagnosticMacAyealPattyn(void); 228 228 ElementMatrix* CreateKMatrixDiagnosticMacAyealStokes(void); 229 ElementMatrix* CreateKMatrixDiagnosticL1L2(void); 230 ElementMatrix* CreateKMatrixDiagnosticL1L2Viscous(void); 231 ElementMatrix* CreateKMatrixDiagnosticL1L2Friction(void); 229 232 ElementMatrix* CreateKMatrixDiagnosticPattyn(void); 230 233 ElementMatrix* CreateKMatrixDiagnosticPattynViscous(void); … … 245 248 void InputUpdateFromSolutionDiagnosticMacAyealPattyn( IssmDouble* solutiong); 246 249 void InputUpdateFromSolutionDiagnosticMacAyealStokes( IssmDouble* solutiong); 250 void InputUpdateFromSolutionDiagnosticL1L2( IssmDouble* solutiong); 247 251 void InputUpdateFromSolutionDiagnosticPattyn( IssmDouble* solutiong); 248 252 void InputUpdateFromSolutionDiagnosticPattynStokes( IssmDouble* solutiong); … … 265 269 ElementVector* CreatePVectorDiagnosticMacAyealPattyn(void); 266 270 ElementVector* CreatePVectorDiagnosticMacAyealStokes(void); 271 ElementVector* CreatePVectorDiagnosticL1L2(void); 267 272 ElementVector* CreatePVectorDiagnosticPattyn(void); 268 273 ElementVector* CreatePVectorDiagnosticPattynStokes(void); … … 273 278 ElementVector* CreatePVectorDiagnosticVertVolume(void); 274 279 ElementVector* CreatePVectorDiagnosticVertBase(void); 280 void GetL1L2Viscosity(IssmDouble*, IssmDouble*, GaussPenta*, Input*, Input*, Input*); 275 281 #endif 276 282 -
issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp
r13036 r13051 232 232 /*Intermediaries */ 233 233 int i,j,ig; 234 IssmDouble 235 IssmDouble 236 IssmDouble 237 IssmDouble 234 IssmDouble heatcapacity,latentheat; 235 IssmDouble Jdet,D_scalar; 236 IssmDouble xyz_list[NUMVERTICES][3]; 237 IssmDouble L[3]; 238 238 GaussTria *gauss=NULL; 239 239 … … 1444 1444 1445 1445 /*Intermediaries*/ 1446 int i,j;1447 int tria_vertex_ids[3];1446 int i,j; 1447 int tria_vertex_ids[3]; 1448 1448 IssmDouble nodeinputs[3]; 1449 1449 IssmDouble cmmininputs[3]; 1450 1450 IssmDouble cmmaxinputs[3]; 1451 bool control_analysis=false;1452 int num_control_type;1451 bool control_analysis = false; 1452 int num_control_type; 1453 1453 IssmDouble yts; 1454 int num_cm_responses;1454 int num_cm_responses; 1455 1455 1456 1456 /*Get parameters: */ -
issm/trunk-jpl/src/c/classes/objects/Materials/Matice.cpp
r13036 r13051 128 128 } 129 129 /*}}}*/ 130 /*FUNCTION Matice::GetA {{{*/ 131 IssmDouble Matice::GetA(){ 132 /* 133 * A = 1/B^n 134 */ 135 136 IssmDouble B,n; 137 138 inputs->GetInputAverage(&B,MaterialsRheologyBEnum); 139 n=this->GetN(); 140 141 return pow(B,-n); 142 } 143 /*}}}*/ 130 144 /*FUNCTION Matice::GetB {{{*/ 131 145 IssmDouble Matice::GetB(){ … … 236 250 if(A==0){ 237 251 /*Maxiviscositym viscosity for 0 shear areas: */ 238 viscosity=2.5*pow( (IssmDouble)10,(IssmDouble)17);252 viscosity=2.5*pow(10.,17.); 239 253 } 240 254 else{ -
issm/trunk-jpl/src/c/classes/objects/Materials/Matice.h
r12472 r13051 58 58 /*}}}*/ 59 59 /*Matice Numerics: {{{*/ 60 void SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin); 61 void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon); 62 void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon); 63 void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon); 64 void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); 65 void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); 66 void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); 60 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); 61 void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon); 62 void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon); 63 void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon); 64 void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); 65 void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); 66 void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); 67 IssmDouble GetA(); 67 68 IssmDouble GetB(); 68 69 IssmDouble GetBbar(); 69 70 IssmDouble GetN(); 70 bool IsInput(int name);71 bool IsInput(int name); 71 72 /*}}}*/ 72 73 }; -
issm/trunk-jpl/src/c/classes/objects/Node.cpp
r13036 r13051 92 92 for(k=1;k<=gsize;k++) this->FreezeDof(k); 93 93 } 94 if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==L1L2ApproximationEnum && !iomodel->Data(MeshVertexonbedEnum)[io_index]){ 95 for(k=1;k<=gsize;k++) this->FreezeDof(k); 96 } 94 97 if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==MacAyealPattynApproximationEnum && iomodel->Data(FlowequationBordermacayealEnum)[io_index]){ 95 98 if(!iomodel->Data(MeshVertexonbedEnum)[io_index]){ … … 114 117 /*Diagnostic Hutter*/ 115 118 if (analysis_type==DiagnosticHutterAnalysisEnum){ 116 if (!iomodel->Data(FlowequationVertexEquationEnum)) _error_("iomodel->vertices_type is NULL");119 _assert_(iomodel->Data(FlowequationVertexEquationEnum)); 117 120 /*Constrain all nodes that are not Hutter*/ 118 121 if (reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[io_index])!=HutterApproximationEnum){ -
issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp
r12832 r13051 36 36 parameters->AddObject(iomodel->CopyConstantObject(FlowequationIshutterEnum)); 37 37 parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsmacayealpattynEnum)); 38 parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsl1l2Enum)); 38 39 parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsstokesEnum)); 39 40 parameters->AddObject(iomodel->CopyConstantObject(SettingsOutputFrequencyEnum)); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp
r13036 r13051 20 20 double rho_ice; 21 21 double stokesreconditioning; 22 bool isstokes,is macayealpattyn;22 bool isstokes,isl1l2,ismacayealpattyn; 23 23 bool spcpresent=false; 24 24 int Mx,Nx; … … 56 56 iomodel->Constant(&stokesreconditioning,DiagnosticStokesreconditioningEnum); 57 57 iomodel->Constant(&isstokes,FlowequationIsstokesEnum); 58 iomodel->Constant(&isl1l2,FlowequationIsl1l2Enum); 58 59 iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum); 59 60 … … 65 66 66 67 /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ 67 if (!ismacayealpattyn & !isstokes){68 if(!ismacayealpattyn & !isstokes & !isl1l2){ 68 69 *pconstraints=constraints; 69 70 return; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp
r12832 r13051 94 94 count++; 95 95 } 96 else if ((int)*(elements_type+element)==(L1L2ApproximationEnum)){ 97 loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,L1L2IceFrontEnum,DiagnosticHorizAnalysisEnum)); 98 count++; 99 } 96 100 else if ((int)*(elements_type+element)==(StokesApproximationEnum)){ 97 101 loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,StokesIceFrontEnum,DiagnosticHorizAnalysisEnum)); -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp
r12832 r13051 19 19 bool continuous_galerkin=true; 20 20 int numberofvertices; 21 bool isstokes,is macayealpattyn;21 bool isstokes,isl1l2,ismacayealpattyn; 22 22 23 23 /*DataSets: */ … … 27 27 iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum); 28 28 iomodel->Constant(&isstokes,FlowequationIsstokesEnum); 29 iomodel->Constant(&isl1l2,FlowequationIsl1l2Enum); 29 30 iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum); 30 31 … … 36 37 37 38 /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ 38 if(!ismacayealpattyn & !isstokes ){39 if(!ismacayealpattyn & !isstokes & !isl1l2){ 39 40 *pnodes=nodes; 40 41 return; -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp
r12832 r13051 19 19 int numberofelements; 20 20 bool ismacayealpattyn; 21 bool isl1l2; 21 22 bool isstokes; 22 23 bool control_analysis; … … 25 26 /*Fetch constants needed: */ 26 27 iomodel->Constant(&isstokes,FlowequationIsstokesEnum); 28 iomodel->Constant(&isl1l2,FlowequationIsl1l2Enum); 27 29 iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum); 28 30 iomodel->Constant(&dim,MeshDimensionEnum); … … 32 34 33 35 /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ 34 if(!ismacayealpattyn & !isstokes ) return;36 if(!ismacayealpattyn & !isstokes &!isl1l2) return; 35 37 36 38 /*Fetch data needed: */ -
issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp
r13036 r13051 19 19 if (analysis_type==DiagnosticHorizAnalysisEnum){ 20 20 if (vertices_type[0]==MacAyealApproximationEnum){ 21 numdofs=2; 22 } 23 else if (vertices_type[0]==L1L2ApproximationEnum){ 21 24 numdofs=2; 22 25 } -
issm/trunk-jpl/src/c/solutions/convergence.cpp
r13036 r13051 81 81 //print 82 82 if(res<eps_res){ 83 if(VerboseConvergence()) _pprintLine_(setw(50) << left << " mechanical equilibrium convergence criterion" << res*100 << " < " << eps_res*100 <<" %");83 if(VerboseConvergence()) _pprintLine_(setw(50)<<left<<" mechanical equilibrium convergence criterion"<<res*100<< " < "<<eps_res*100<<" %"); 84 84 converged=true; 85 85 } 86 86 else{ 87 if(VerboseConvergence()) _pprintLine_(setw(50) << left << " mechanical equilibrium convergence criterion" << res*100 << " > " << eps_res*100 <<" %");87 if(VerboseConvergence()) _pprintLine_(setw(50)<<left<<" mechanical equilibrium convergence criterion"<<res*100<<" > "<<eps_res*100<<" %"); 88 88 converged=false; 89 89 } -
issm/trunk-jpl/src/c/solutions/diagnostic_core.cpp
r12832 r13051 20 20 bool ishutter = false; 21 21 bool ismacayealpattyn = false; 22 bool isl1l2 = false; 22 23 bool isstokes = false; 23 24 bool isnewton = false; … … 33 34 femmodel->parameters->FindParam(&ishutter,FlowequationIshutterEnum); 34 35 femmodel->parameters->FindParam(&ismacayealpattyn,FlowequationIsmacayealpattynEnum); 36 femmodel->parameters->FindParam(&isl1l2,FlowequationIsl1l2Enum); 35 37 femmodel->parameters->FindParam(&isstokes,FlowequationIsstokesEnum); 36 38 femmodel->parameters->FindParam(&isnewton,DiagnosticIsnewtonEnum); … … 70 72 } 71 73 72 if ( ismacayealpattyn^ isstokes){ // ^ = xor74 if ((ismacayealpattyn || isl1l2) ^ isstokes){ // ^ = xor 73 75 74 76 if(VerboseSolution()) _pprintLine_(" computing velocities");
Note:
See TracChangeset
for help on using the changeset viewer.