Ignore:
Timestamp:
11/15/13 10:49:47 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: InputUpdateFromSolution is now entirely done by analyses

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r16775 r16783  
    22742274                this->inputs->AddInput(datasetinput);
    22752275        }
    2276 }
    2277 /*}}}*/
    2278 /*FUNCTION Penta::InputUpdateFromSolution {{{*/
    2279 void  Penta::InputUpdateFromSolution(IssmDouble* solution){
    2280 
    2281         int analysis_type,inputenum;
    2282 
    2283         /*retreive parameters: */
    2284         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    2285 
    2286         /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
    2287         switch(analysis_type){
    2288         #ifdef _HAVE_STRESSBALANCE_
    2289         case StressbalanceAnalysisEnum:
    2290                 InputUpdateFromSolutionStressbalanceHoriz( solution);
    2291                 break;
    2292         case StressbalanceSIAAnalysisEnum:
    2293                 InputUpdateFromSolutionStressbalanceSIA( solution);
    2294                 break;
    2295         case StressbalanceVerticalAnalysisEnum:
    2296                 InputUpdateFromSolutionStressbalanceVert( solution);
    2297                 break;
    2298         #endif
    2299         #ifdef _HAVE_CONTROL_
    2300         case AdjointHorizAnalysisEnum:
    2301                 int approximation;
    2302                 inputs->GetInputValue(&approximation,ApproximationEnum);
    2303                 if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
    2304                         InputUpdateFromSolutionAdjointFS( solution);
    2305                 }
    2306                 else{
    2307                         InputUpdateFromSolutionAdjointHoriz( solution);
    2308                 }
    2309                 break;
    2310         #endif
    2311         #ifdef _HAVE_THERMAL_
    2312         case ThermalAnalysisEnum:
    2313                 InputUpdateFromSolutionThermal( solution);
    2314                 break;
    2315         case EnthalpyAnalysisEnum:
    2316                 InputUpdateFromSolutionEnthalpy( solution);
    2317                 break;
    2318         case MeltingAnalysisEnum:
    2319                 InputUpdateFromSolutionOneDof(solution,BasalforcingsMeltingRateEnum);
    2320                 break;
    2321         #endif
    2322         case L2ProjectionBaseAnalysisEnum:
    2323                 this->parameters->FindParam(&inputenum,InputToL2ProjectEnum);
    2324                 InputUpdateFromSolutionOneDofCollapsed(solution,inputenum);
    2325                 break;
    2326         case MasstransportAnalysisEnum:
    2327                 InputUpdateFromSolutionMasstransport(solution);
    2328                 break;
    2329         case FreeSurfaceTopAnalysisEnum:
    2330                 InputUpdateFromSolutionFreeSurfaceTop(solution);
    2331                 break;
    2332         case FreeSurfaceBaseAnalysisEnum:
    2333                 InputUpdateFromSolutionFreeSurfaceBase(solution);
    2334                 break;
    2335         #ifdef _HAVE_BALANCED_
    2336         case BalancethicknessAnalysisEnum:
    2337                 InputUpdateFromSolutionOneDofCollapsed(solution,ThicknessEnum);
    2338                 break;
    2339         #endif
    2340         #ifdef _HAVE_HYDROLOGY_
    2341         case HydrologyDCInefficientAnalysisEnum:
    2342                 InputUpdateFromSolutionHydrologyDCInefficient(solution);
    2343                 break;
    2344         case HydrologyDCEfficientAnalysisEnum:
    2345                 InputUpdateFromSolutionOneDofCollapsed(solution,EplHeadEnum);
    2346                 break;
    2347         case L2ProjectionEPLAnalysisEnum:
    2348                 this->parameters->FindParam(&inputenum,InputToL2ProjectEnum);
    2349                 InputUpdateFromSolutionOneDofCollapsed(solution,inputenum);
    2350                 break;
    2351         #endif
    2352         default:
    2353                 _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    2354         }
    2355 }
    2356 /*}}}*/
    2357 /*FUNCTION Penta::InputUpdateFromSolutionMasstransport{{{*/
    2358 void  Penta::InputUpdateFromSolutionMasstransport(IssmDouble* solution){
    2359 
    2360         const int  numdof   = NDOF1*NUMVERTICES;
    2361         const int  numdof2d = NDOF1*NUMVERTICES2D;
    2362 
    2363         int    i,hydroadjustment;
    2364         int*   doflist = NULL;
    2365         IssmDouble rho_ice,rho_water,minthickness;
    2366         IssmDouble newthickness[numdof];
    2367         IssmDouble newbed[numdof];
    2368         IssmDouble newsurface[numdof];
    2369         IssmDouble oldbed[NUMVERTICES];
    2370         IssmDouble oldsurface[NUMVERTICES];
    2371         IssmDouble oldthickness[NUMVERTICES];
    2372         IssmDouble phi[NUMVERTICES];
    2373         Penta  *penta   = NULL;
    2374 
    2375         /*If not on bed, return*/
    2376         if (!IsOnBed()) return;
    2377 
    2378         /*Get dof list: */
    2379         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    2380 
    2381         /*Use the dof list to index into the solution vector and extrude it */
    2382         this->parameters->FindParam(&minthickness,MasstransportMinThicknessEnum);
    2383         for(i=0;i<numdof2d;i++){
    2384                 newthickness[i]=solution[doflist[i]];
    2385                 if(xIsNan<IssmDouble>(newthickness[i])) _error_("NaN found in solution vector");
    2386                 /*Constrain thickness to be at least 1m*/
    2387                 if(newthickness[i]<minthickness) newthickness[i]=minthickness;
    2388                 newthickness[i+numdof2d]=newthickness[i];
    2389         }
    2390 
    2391         /*Get previous bed, thickness and surface*/
    2392         GetInputListOnVertices(&oldbed[0],BedEnum);
    2393         GetInputListOnVertices(&oldsurface[0],SurfaceEnum);
    2394         GetInputListOnVertices(&oldthickness[0],ThicknessEnum);
    2395         GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
    2396 
    2397         /*Fing MasstransportHydrostaticAdjustment to figure out how to update the geometry:*/
    2398         this->parameters->FindParam(&hydroadjustment,MasstransportHydrostaticAdjustmentEnum);
    2399 
    2400         /*recover material parameters: */
    2401         rho_ice=matpar->GetRhoIce();
    2402         rho_water=matpar->GetRhoWater();
    2403 
    2404         for(i=0;i<numdof;i++) {
    2405                 /*If shelf: hydrostatic equilibrium*/
    2406                 if (phi[i]>0.){
    2407                         newsurface[i]=oldbed[i]+newthickness[i]; //surface = oldbed + newthickness
    2408                         newbed[i]=oldbed[i];               //same bed: do nothing
    2409                 }
    2410                 else{ //so it is an ice shelf
    2411                         if(hydroadjustment==AbsoluteEnum){
    2412                                 newsurface[i]=newthickness[i]*(1-rho_ice/rho_water);
    2413                                 newbed[i]=newthickness[i]*(-rho_ice/rho_water);
    2414                         }
    2415                         else if(hydroadjustment==IncrementalEnum){
    2416                                 newsurface[i]=oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH
    2417                                 newbed[i]=oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed = oldbed + di * dH
    2418                         }
    2419                         else _error_("Hydrostatic adjustment " << hydroadjustment << " (" << EnumToStringx(hydroadjustment) << ") not supported yet");
    2420                 }
    2421         }
    2422 
    2423         /*Start looping over all elements above current element and update all inputs*/
    2424         penta=this;
    2425         for(;;){
    2426                 /*Add input to the element: */
    2427                 penta->inputs->AddInput(new PentaInput(ThicknessEnum,newthickness,P1Enum));
    2428                 penta->inputs->AddInput(new PentaInput(SurfaceEnum,newsurface,P1Enum));
    2429                 penta->inputs->AddInput(new PentaInput(BedEnum,newbed,P1Enum));
    2430 
    2431                 /*Stop if we have reached the surface*/
    2432                 if (penta->IsOnSurface()) break;
    2433 
    2434                 /* get upper Penta*/
    2435                 penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
    2436         }
    2437 
    2438         /*Free ressources:*/
    2439         xDelete<int>(doflist);
    2440 }
    2441 /*}}}*/
    2442 /*FUNCTION Penta::InputUpdateFromSolutionFreeSurfaceTop{{{*/
    2443 void  Penta::InputUpdateFromSolutionFreeSurfaceTop(IssmDouble* solution){
    2444 
    2445         const int  numdof   = NDOF1*NUMVERTICES;
    2446         const int  numdof2d = NDOF1*NUMVERTICES2D;
    2447 
    2448         int    i;
    2449         int*   doflist = NULL;
    2450         IssmDouble newthickness[numdof];
    2451         IssmDouble newbed[numdof];
    2452         IssmDouble newsurface[numdof];
    2453 
    2454         /*If not on bed, return*/
    2455         if (!IsOnSurface()) return;
    2456 
    2457         /*Get dof list: */
    2458         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    2459 
    2460         /*Use the dof list to index into the solution vector and extrude it */
    2461         for(i=0;i<numdof2d;i++){
    2462                 newsurface[i+numdof2d]=solution[doflist[i+numdof2d]];
    2463                 if(xIsNan<IssmDouble>(newsurface[i+numdof2d])) _error_("NaN found in solution vector");
    2464                 newsurface[i]=newsurface[i+numdof2d];
    2465         }
    2466 
    2467         /*Get previous bed and thickness*/
    2468         GetInputListOnVertices(&newbed[0],BedEnum);
    2469 
    2470         for(i=0;i<numdof;i++) {
    2471                 newthickness[i]=newsurface[i]-newbed[i];
    2472                 _assert_(newthickness[i]>0.);
    2473         }
    2474 
    2475         /*Start looping over all elements above current element and update all inputs*/
    2476         Penta* penta=this;
    2477         for(;;){
    2478                 /*Add input to the element: */
    2479                 penta->inputs->AddInput(new PentaInput(ThicknessEnum,newthickness,P1Enum));
    2480                 penta->inputs->AddInput(new PentaInput(SurfaceEnum,newsurface,P1Enum));
    2481 
    2482                 /*Stop if we have reached the surface*/
    2483                 if(penta->IsOnBed()) break;
    2484 
    2485                 /* get upper Penta*/
    2486                 penta=penta->GetLowerElement(); _assert_(penta->Id()!=this->id);
    2487         }
    2488 
    2489         /*Free ressources:*/
    2490         xDelete<int>(doflist);
    2491 }
    2492 /*}}}*/
    2493 /*FUNCTION Penta::InputUpdateFromSolutionFreeSurfaceBase{{{*/
    2494 void  Penta::InputUpdateFromSolutionFreeSurfaceBase(IssmDouble* solution){
    2495 
    2496         const int  numdof   = NDOF1*NUMVERTICES;
    2497         const int  numdof2d = NDOF1*NUMVERTICES2D;
    2498 
    2499         int    i;
    2500         int*   doflist = NULL;
    2501         IssmDouble newthickness[numdof];
    2502         IssmDouble newbed[numdof];
    2503         IssmDouble newsurface[numdof];
    2504 
    2505         /*If not on bed, return*/
    2506         if (!IsOnBed()) return;
    2507 
    2508         /*Get dof list: */
    2509         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    2510 
    2511         /*Use the dof list to index into the solution vector and extrude it */
    2512         for(i=0;i<numdof2d;i++){
    2513                 newbed[i]=solution[doflist[i]];
    2514                 if(xIsNan<IssmDouble>(newbed[i])) _error_("NaN found in solution vector");
    2515                 newbed[i+numdof2d]=newbed[i];
    2516         }
    2517 
    2518         /*Get previous bed and thickness*/
    2519         GetInputListOnVertices(&newsurface[0],SurfaceEnum);
    2520 
    2521         for(i=0;i<numdof;i++) {
    2522                 newthickness[i]=newsurface[i]-newbed[i];
    2523                 _assert_(newthickness[i]>0.);
    2524         }
    2525 
    2526         /*Start looping over all elements above current element and update all inputs*/
    2527         Penta* penta=this;
    2528         for(;;){
    2529                 /*Add input to the element: */
    2530                 penta->inputs->AddInput(new PentaInput(ThicknessEnum,newthickness,P1Enum));
    2531                 penta->inputs->AddInput(new PentaInput(BedEnum,newbed,P1Enum));
    2532 
    2533                 /*Stop if we have reached the surface*/
    2534                 if(penta->IsOnSurface()) break;
    2535 
    2536                 /* get upper Penta*/
    2537                 penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
    2538         }
    2539 
    2540         /*Free ressources:*/
    2541         xDelete<int>(doflist);
    25422276}
    25432277/*}}}*/
    … …  
    48544588        delete friction;
    48554589        return pe;
    4856 }
    4857 /*}}}*/
    4858 /*FUNCTION Penta::InputUpdateFromSolutionThermal {{{*/
    4859 void  Penta::InputUpdateFromSolutionThermal(IssmDouble* solution){
    4860 
    4861         const int    numdof=NDOF1*NUMVERTICES;
    4862 
    4863         bool        converged;
    4864         int         i,rheology_law;
    4865         IssmDouble  xyz_list[NUMVERTICES][3];
    4866         IssmDouble  values[numdof];
    4867         IssmDouble  B[numdof],surface[numdof];
    4868         IssmDouble  B_average;
    4869         int        *doflist = NULL;
    4870         bool        hack    = false;
    4871 
    4872         /*Get dof list: */
    4873         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    4874 
    4875         /*Use the dof list to index into the solution vector: */
    4876         for(i=0;i<numdof;i++){
    4877                 values[i]=solution[doflist[i]];
    4878 
    4879                 /*Check solution*/
    4880                 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
    4881                 //if(values[i]<0)      _printf_("temperature < 0°K found in solution vector\n");
    4882                 //if(values[i]>275)    _printf_("temperature > 275°K found in solution vector (Paterson's rheology associated is negative)\n");
    4883         }
    4884 
    4885         if(hack){
    4886                 /*Force temperature between [Tpmp-50 Tpmp] to disable penalties*/
    4887                 IssmDouble pressure[numdof];
    4888                 GetInputListOnVertices(&pressure[0],PressureEnum);
    4889                 for(i=0;i<numdof;i++){
    4890                         if(values[i]>matpar->TMeltingPoint(pressure[i]))     values[i]=matpar->TMeltingPoint(pressure[i]);
    4891                         if(values[i]<matpar->TMeltingPoint(pressure[i])-50.) values[i]=matpar->TMeltingPoint(pressure[i])-50.;
    4892                 }
    4893         }
    4894 
    4895 
    4896         /*Get all inputs and parameters*/
    4897         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4898         Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
    4899 
    4900         this->inputs->GetInputValue(&converged,ConvergedEnum);
    4901         if(converged){
    4902                 this->inputs->AddInput(new PentaInput(TemperatureEnum,values,P1Enum));
    4903 
    4904                 /*Update Rheology only if converged (we must make sure that the temperature is below melting point
    4905                  * otherwise the rheology could be negative*/
    4906                 this->parameters->FindParam(&rheology_law,MaterialsRheologyLawEnum);
    4907                 GetInputListOnVertices(&surface[0],SurfaceEnum);
    4908                 switch(rheology_law){
    4909                         case NoneEnum:
    4910                                 /*Do nothing: B is not temperature dependent*/
    4911                                 break;
    4912                         case PatersonEnum:
    4913                                 for(i=0;i<numdof;i++) B[i]=Paterson(values[i]);
    4914                                 this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
    4915                                 break;
    4916                         case ArrheniusEnum:
    4917                                 for(i=0;i<numdof;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i][2],material->GetN());
    4918                                 this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
    4919                                 break;
    4920                         default:
    4921                                 _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
    4922 
    4923                 }
    4924         }
    4925         else{
    4926                 this->inputs->AddInput(new PentaInput(TemperaturePicardEnum,values,P1Enum));
    4927         }
    4928 
    4929         /*Free ressources:*/
    4930         xDelete<int>(doflist);
    4931 }
    4932 /*}}}*/
    4933 /*FUNCTION Penta::InputUpdateFromSolutionEnthalpy {{{*/
    4934 void  Penta::InputUpdateFromSolutionEnthalpy(IssmDouble* solution){
    4935 
    4936         const int    numdof=NDOF1*NUMVERTICES;
    4937 
    4938         bool       converged=false;
    4939         int        i,rheology_law;
    4940         IssmDouble xyz_list[NUMVERTICES][3];
    4941         IssmDouble values[numdof];
    4942         IssmDouble pressure[numdof];
    4943         IssmDouble surface[numdof];
    4944         IssmDouble temperatures[numdof];
    4945         IssmDouble waterfraction[numdof];
    4946         IssmDouble B[numdof];
    4947         int*       doflist=NULL;
    4948 
    4949         /*Get dof list: */
    4950         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    4951 
    4952         /*Use the dof list to index into the solution vector: */
    4953         for(i=0;i<numdof;i++){
    4954                 values[i]=solution[doflist[i]];
    4955 
    4956                 /*Check solution*/
    4957                 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
    4958         }
    4959 
    4960         /*Get all inputs and parameters*/
    4961         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4962         GetInputListOnVertices(&pressure[0],PressureEnum);
    4963         Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
    4964 
    4965         this->inputs->GetInputValue(&converged,ConvergedEnum);
    4966         if(converged){
    4967                 /*Convert enthalpy into temperature and water fraction*/
    4968                 for(i=0;i<numdof;i++){
    4969                         matpar->EnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]);
    4970                         if(waterfraction[i]<0) _error_("Negative water fraction found in solution vector");
    4971                         //if(waterfraction[i]>1) _error_("Water fraction >1 found in solution vector");
    4972                 }
    4973 
    4974                 this->inputs->AddInput(new PentaInput(EnthalpyEnum,values,P1Enum));
    4975                 this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum));
    4976                 this->inputs->AddInput(new PentaInput(TemperatureEnum,temperatures,P1Enum));
    4977 
    4978                 /*Update Rheology only if converged (we must make sure that the temperature is below melting point
    4979                  * otherwise the rheology could be negative*/
    4980                 this->parameters->FindParam(&rheology_law,MaterialsRheologyLawEnum);
    4981                 GetInputListOnVertices(&surface[0],SurfaceEnum);
    4982                 switch(rheology_law){
    4983                         case NoneEnum:
    4984                                 /*Do nothing: B is not temperature dependent*/
    4985                                 break;
    4986                         case PatersonEnum:
    4987                                 for(i=0;i<numdof;i++) B[i]=Paterson(temperatures[i]);
    4988                                 this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
    4989                                 break;
    4990                         case ArrheniusEnum:
    4991                                 for(i=0;i<numdof;i++) B[i]=Arrhenius(temperatures[i],surface[i]-xyz_list[i][2],material->GetN());
    4992                                 this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
    4993                                 break;
    4994                         case LliboutryDuvalEnum:
    4995                                 for(i=0;i<numdof;i++) B[i]=LliboutryDuval(values[i],pressure[i],material->GetN(),matpar->GetBeta(),matpar->GetReferenceTemperature(),matpar->GetHeatCapacity(),matpar->GetLatentHeat());
    4996                                 this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
    4997                                 break;
    4998                         default:
    4999                                 _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
    5000                 }
    5001         }
    5002         else{
    5003                 this->inputs->AddInput(new PentaInput(EnthalpyPicardEnum,values,P1Enum));
    5004         }
    5005 
    5006         /*Free ressources:*/
    5007         xDelete<int>(doflist);
    50084590}
    50094591/*}}}*/
    … …  
    63715953}
    63725954/*}}}*/
    6373 /*FUNCTION Penta::InputUpdateFromSolutionAdjointFS {{{*/
    6374 void  Penta::InputUpdateFromSolutionAdjointFS(IssmDouble* solution){
    6375 
    6376         int          i;
    6377         int*         vdoflist=NULL;
    6378         int*         pdoflist=NULL;
    6379         IssmDouble   FSreconditioning;
    6380 
    6381         /*Fetch number of nodes and dof for this finite element*/
    6382         int vnumnodes = this->NumberofNodesVelocity();
    6383         int pnumnodes = this->NumberofNodesPressure();
    6384         int vnumdof   = vnumnodes*NDOF3;
    6385         int pnumdof   = pnumnodes*NDOF1;
    6386 
    6387         /*Initialize values*/
    6388         IssmDouble* vvalues = xNew<IssmDouble>(vnumdof);
    6389         IssmDouble* pvalues = xNew<IssmDouble>(pnumdof);
    6390         IssmDouble* lambdax = xNew<IssmDouble>(vnumnodes);
    6391         IssmDouble* lambday = xNew<IssmDouble>(vnumnodes);
    6392         IssmDouble* lambdaz = xNew<IssmDouble>(vnumnodes);
    6393         IssmDouble* lambdap = xNew<IssmDouble>(pnumnodes);
    6394 
    6395         /*Get dof list: */
    6396         GetDofListVelocity(&vdoflist,GsetEnum);
    6397         GetDofListPressure(&pdoflist,GsetEnum);
    6398 
    6399         /*Use the dof list to index into the solution vector: */
    6400         for(i=0;i<vnumdof;i++) vvalues[i]=solution[vdoflist[i]];
    6401         for(i=0;i<pnumdof;i++) pvalues[i]=solution[pdoflist[i]];
    6402 
    6403         /*Transform solution in Cartesian Space*/
    6404         ::TransformSolutionCoord(&vvalues[0],nodes,vnumnodes,XYZEnum);
    6405 
    6406         /*fill in all arrays: */
    6407         for(i=0;i<vnumnodes;i++){
    6408                 lambdax[i] = vvalues[i*NDOF3+0]; if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector");
    6409                 lambday[i] = vvalues[i*NDOF3+1]; if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
    6410                 lambdaz[i] = vvalues[i*NDOF3+2]; if(xIsNan<IssmDouble>(lambdaz[i])) _error_("NaN found in solution vector");
    6411         }
    6412         for(i=0;i<pnumnodes;i++){
    6413                 lambdap[i] = pvalues[i]; if(xIsNan<IssmDouble>(lambdap[i])) _error_("NaN found in solution vector");
    6414         }
    6415 
    6416         /*Recondition pressure and compute vel: */
    6417         this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    6418         for(i=0;i<pnumnodes;i++) lambdap[i]=lambdap[i]*FSreconditioning;
    6419 
    6420         /*Add vx and vy as inputs to the tria element: */
    6421         this->inputs->AddInput(new PentaInput(AdjointxEnum,lambdax,P1Enum));
    6422         this->inputs->AddInput(new PentaInput(AdjointyEnum,lambday,P1Enum));
    6423         this->inputs->AddInput(new PentaInput(AdjointzEnum,lambdaz,P1Enum));
    6424         this->inputs->AddInput(new PentaInput(AdjointpEnum,lambdap,P1Enum));
    6425 
    6426         /*Free ressources:*/
    6427         xDelete<int>(vdoflist);
    6428         xDelete<int>(pdoflist);
    6429         xDelete<IssmDouble>(lambdap);
    6430         xDelete<IssmDouble>(lambdaz);
    6431         xDelete<IssmDouble>(lambday);
    6432         xDelete<IssmDouble>(lambdax);
    6433         xDelete<IssmDouble>(vvalues);
    6434         xDelete<IssmDouble>(pvalues);
    6435 }
    6436 /*}}}*/
    6437 /*FUNCTION Penta::InputUpdateFromSolutionAdjointHoriz {{{*/
    6438 void  Penta::InputUpdateFromSolutionAdjointHoriz(IssmDouble* solution){
    6439 
    6440         const int numdof=NDOF2*NUMVERTICES;
    6441 
    6442         int    i;
    6443         IssmDouble values[numdof];
    6444         IssmDouble lambdax[NUMVERTICES];
    6445         IssmDouble lambday[NUMVERTICES];
    6446         int*   doflist=NULL;
    6447 
    6448         /*Get dof list: */
    6449         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    6450 
    6451         /*Use the dof list to index into the solution vector: */
    6452         for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    6453 
    6454         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    6455         for(i=0;i<NUMVERTICES;i++){
    6456                 lambdax[i]=values[i*NDOF2+0];
    6457                 lambday[i]=values[i*NDOF2+1];
    6458 
    6459                 /*Check solution*/
    6460                 if(xIsNan<IssmDouble>(lambdax[i]))       _error_("NaN found in solution vector");
    6461                 if(xIsNan<IssmDouble>(lambday[i]))       _error_("NaN found in solution vector");
    6462         }
    6463 
    6464         /*Add vx and vy as inputs to the tria element: */
    6465         this->inputs->AddInput(new PentaInput(AdjointxEnum,lambdax,P1Enum));
    6466         this->inputs->AddInput(new PentaInput(AdjointyEnum,lambday,P1Enum));
    6467 
    6468         /*Free ressources:*/
    6469         xDelete<int>(doflist);
    6470 }
    6471 /*}}}*/
    64725955/*FUNCTION Penta::SurfaceAverageVelMisfit {{{*/
    64735956IssmDouble Penta::SurfaceAverageVelMisfit(void){
    … …  
    98469329        *pviscosity = viscosity;
    98479330        return;
    9848 }
    9849 /*}}}*/
    9850 /*FUNCTION Penta::InputUpdateFromSolutionStressbalanceHoriz {{{*/
    9851 void  Penta::InputUpdateFromSolutionStressbalanceHoriz(IssmDouble* solution){
    9852 
    9853         int  approximation;
    9854 
    9855         /*Recover inputs*/
    9856         inputs->GetInputValue(&approximation,ApproximationEnum);
    9857 
    9858         /*SSA, everything is done by the element on bed*/
    9859         if (approximation==SSAApproximationEnum){
    9860                 if (!IsOnBed()){
    9861                         /*Do nothing. Element on bed will take care of it*/
    9862                         return;
    9863                 }
    9864                 else{
    9865                         InputUpdateFromSolutionStressbalanceSSA(solution);
    9866                         return;
    9867                 }
    9868         }
    9869         if (approximation==L1L2ApproximationEnum){
    9870                 if (!IsOnBed()) return;
    9871                 InputUpdateFromSolutionStressbalanceL1L2(solution);
    9872                 return;
    9873         }
    9874         else if (approximation==HOApproximationEnum){
    9875                 InputUpdateFromSolutionStressbalanceHO(solution);
    9876         }
    9877         else if (approximation==HOFSApproximationEnum){
    9878                 InputUpdateFromSolutionStressbalanceHOFS(solution);
    9879         }
    9880         else if (approximation==SSAFSApproximationEnum){
    9881                 InputUpdateFromSolutionStressbalanceSSAFS(solution);
    9882         }
    9883         else if (approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
    9884                 InputUpdateFromSolutionStressbalanceFS(solution);
    9885         }
    9886         else if (approximation==SSAHOApproximationEnum){
    9887                 InputUpdateFromSolutionStressbalanceSSAHO(solution);
    9888         }
    9889 }
    9890 /*}}}*/
    9891 /*FUNCTION Penta::InputUpdateFromSolutionStressbalanceSSA {{{*/
    9892 void  Penta::InputUpdateFromSolutionStressbalanceSSA(IssmDouble* solution){
    9893 
    9894         int         numnodes = this->NumberofNodes();
    9895         int         numdof   = NDOF2*numnodes;
    9896 
    9897         int         i;
    9898         IssmDouble  rho_ice,g;
    9899         IssmDouble  values[2*NUMVERTICES];
    9900         IssmDouble  vx[NUMVERTICES];
    9901         IssmDouble  vy[NUMVERTICES];
    9902         IssmDouble  vz[NUMVERTICES];
    9903         IssmDouble  vel[NUMVERTICES];
    9904         IssmDouble  pressure[NUMVERTICES];
    9905         IssmDouble  surface[NUMVERTICES];
    9906         IssmDouble  xyz_list[NUMVERTICES][3];
    9907         int        *doflist = NULL;
    9908         Penta      *penta   = NULL;
    9909 
    9910         /*Get dof list: */
    9911         GetDofList(&doflist,SSAApproximationEnum,GsetEnum);
    9912 
    9913         /*Use the dof list to index into the solution vector: */
    9914         for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    9915 
    9916         /*Transform solution in Cartesian Space*/
    9917         ::TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,XYEnum); /*2D: only the first 3 nodes are taken*/
    9918 
    9919         /*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */
    9920         for(i=0;i<3;i++){
    9921                 vx[i]  =values[i*NDOF2+0];
    9922                 vy[i]  =values[i*NDOF2+1];
    9923                 vx[i+3]=vx[i];
    9924                 vy[i+3]=vy[i];
    9925 
    9926                 /*Check solution*/
    9927                 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    9928                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
    9929         }
    9930 
    9931         /*Get parameters fro pressure computation*/
    9932         rho_ice=matpar->GetRhoIce();
    9933         g=matpar->GetG();
    9934 
    9935         /*Start looping over all elements above current element and update all inputs*/
    9936         penta=this;
    9937         for(;;){
    9938 
    9939                 /*Get node data: */
    9940                 ::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
    9941 
    9942                 /*Now Compute vel*/
    9943                 GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
    9944                 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);
    9945 
    9946                 /*Now compute pressure*/
    9947                 GetInputListOnVertices(&surface[0],SurfaceEnum);
    9948                 for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
    9949 
    9950                 /*Now, we have to move the previous Vx and Vy inputs  to old
    9951                  * status, otherwise, we'll wipe them off: */
    9952                 penta->inputs->ChangeEnum(VxEnum,VxPicardEnum);
    9953                 penta->inputs->ChangeEnum(VyEnum,VyPicardEnum);
    9954                 penta->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
    9955 
    9956                 /*Add vx and vy as inputs to the tria element: */
    9957                 penta->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
    9958                 penta->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
    9959                 penta->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
    9960                 penta->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    9961 
    9962                 /*Stop if we have reached the surface*/
    9963                 if (penta->IsOnSurface()) break;
    9964 
    9965                 /* get upper Penta*/
    9966                 penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
    9967         }
    9968 
    9969         /*Free ressources:*/
    9970         xDelete<int>(doflist);
    9971 }
    9972 /*}}}*/
    9973 /*FUNCTION Penta::InputUpdateFromSolutionStressbalanceSSAHO {{{*/
    9974 void  Penta::InputUpdateFromSolutionStressbalanceSSAHO(IssmDouble* solution){
    9975 
    9976         const int    numdof=NDOF2*NUMVERTICES;
    9977         const int    numdof2d=NDOF2*NUMVERTICES2D;
    9978 
    9979         int     i;
    9980         IssmDouble  rho_ice,g;
    9981         IssmDouble  SSA_values[numdof];
    9982         IssmDouble  HO_values[numdof];
    9983         IssmDouble  vx[NUMVERTICES];
    9984         IssmDouble  vy[NUMVERTICES];
    9985         IssmDouble  vz[NUMVERTICES];
    9986         IssmDouble  vel[NUMVERTICES];
    9987         IssmDouble  pressure[NUMVERTICES];
    9988         IssmDouble  surface[NUMVERTICES];
    9989         IssmDouble  xyz_list[NUMVERTICES][3];
    9990         int*    doflistp = NULL;
    9991         int*    doflistm = NULL;
    9992         Penta   *penta   = NULL;
    9993 
    9994         /*OK, we have to add results of this element for HO
    9995          * and results from the penta at base for SSA. Now recover results*/
    9996         penta=(Penta*)GetBasalElement();
    9997 
    9998         /*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
    9999         GetDofList(&doflistp,HOApproximationEnum,GsetEnum);
    10000         penta->GetDofList(&doflistm,SSAApproximationEnum,GsetEnum);
    10001 
    10002         /*Get node data: */
    10003         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    10004 
    10005         /*Use the dof list to index into the solution vector: */
    10006         for(i=0;i<numdof2d;i++){
    10007                 HO_values[i]=solution[doflistp[i]];
    10008                 SSA_values[i]=solution[doflistm[i]];
    10009         }
    10010         for(i=numdof2d;i<numdof;i++){
    10011                 HO_values[i]=solution[doflistp[i]];
    10012                 SSA_values[i]=SSA_values[i-numdof2d];
    10013         }
    10014 
    10015         /*Transform solution in Cartesian Space*/
    10016         ::TransformSolutionCoord(&SSA_values[0],penta->nodes,NUMVERTICES,XYEnum);
    10017         ::TransformSolutionCoord(&HO_values[0],   this->nodes,NUMVERTICES,XYEnum);
    10018 
    10019         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    10020         for(i=0;i<NUMVERTICES;i++){
    10021                 vx[i]=SSA_values[i*NDOF2+0]+HO_values[i*NDOF2+0];
    10022                 vy[i]=SSA_values[i*NDOF2+1]+HO_values[i*NDOF2+1];
    10023 
    10024                 /*Check solution*/
    10025                 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    10026                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
    10027         }
    10028 
    10029         /*Now Compute vel*/
    10030         GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
    10031         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);
    10032 
    10033         /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D,
    10034          *so the pressure is just the pressure at the z elevation: */
    10035         rho_ice=matpar->GetRhoIce();
    10036         g=matpar->GetG();
    10037         GetInputListOnVertices(&surface[0],SurfaceEnum);
    10038         for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
    10039 
    10040         /*Now, we have to move the previous Vx and Vy inputs  to old
    10041          * status, otherwise, we'll wipe them off: */
    10042         this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
    10043         this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
    10044         this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
    10045 
    10046         /*Add vx and vy as inputs to the tria element: */
    10047         this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
    10048         this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
    10049         this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
    10050         this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    10051 
    10052         /*Free ressources:*/
    10053         xDelete<int>(doflistp);
    10054         xDelete<int>(doflistm);
    10055 }
    10056 /*}}}*/
    10057 /*FUNCTION Penta::InputUpdateFromSolutionStressbalanceSSAFS {{{*/
    10058 void  Penta::InputUpdateFromSolutionStressbalanceSSAFS(IssmDouble* solution){
    10059 
    10060         const int    numdofSSA = NDOF2*NUMVERTICES;
    10061         const int    numdof2d  = NDOF2*NUMVERTICES2D;
    10062         const int    numdofFSv = NDOF3*NUMVERTICES;
    10063         const int    numdofFSp = NDOF1*NUMVERTICES;
    10064 
    10065         int     i;
    10066         IssmDouble  FSreconditioning;
    10067         IssmDouble  SSA_values[numdofSSA];
    10068         IssmDouble  FS_values[numdofFSv+numdofFSp];
    10069         IssmDouble  vx[NUMVERTICES];
    10070         IssmDouble  vy[NUMVERTICES];
    10071         IssmDouble  vz[NUMVERTICES];
    10072         IssmDouble  vzSSA[NUMVERTICES];
    10073         IssmDouble  vzFS[NUMVERTICES];
    10074         IssmDouble  vel[NUMVERTICES];
    10075         IssmDouble  pressure[NUMVERTICES];
    10076         IssmDouble  xyz_list[NUMVERTICES][3];
    10077         int   *doflistSSA = NULL;
    10078         int   *doflistFSv = NULL;
    10079         int   *doflistFSp = NULL;
    10080         Penta *penta      = NULL;
    10081 
    10082         /*Prepare coordinate system list*/
    10083         int* cs_list = xNew<int>(2*NUMVERTICES);
    10084         for(i=0;i<NUMVERTICES;i++) cs_list[i]             = XYZEnum;
    10085         for(i=0;i<NUMVERTICES;i++) cs_list[NUMVERTICES+i] = PressureEnum;
    10086 
    10087         /*OK, we have to add results of this element for SSA
    10088          * and results from the penta at base for SSA. Now recover results*/
    10089         penta=(Penta*)GetBasalElement();
    10090 
    10091         /*Get dof listof this element (SSA dofs) and of the penta at base (SSA dofs): */
    10092         penta->GetDofList(&doflistSSA,SSAApproximationEnum,GsetEnum);
    10093         GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
    10094         GetDofListPressure(&doflistFSp,GsetEnum);
    10095         this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    10096 
    10097         /*Get node data: */
    10098         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    10099 
    10100         /*Use the dof list to index into the solution vector: */
    10101         for(i=0;i<numdof2d;i++){
    10102                 SSA_values[i]=solution[doflistSSA[i]];
    10103                 SSA_values[i+numdof2d]=solution[doflistSSA[i]];
    10104         }
    10105         for(i=0;i<numdofFSv;i++) FS_values[i]=solution[doflistFSv[i]];
    10106         for(i=0;i<numdofFSp;i++) FS_values[numdofFSv+i]=solution[doflistFSp[i]];
    10107 
    10108         /*Transform solution in Cartesian Space*/
    10109         ::TransformSolutionCoord(&SSA_values[0],this->nodes,NUMVERTICES,XYEnum);
    10110         ::TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
    10111 
    10112         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    10113         for(i=0;i<NUMVERTICES;i++){
    10114                 vx[i]       = FS_values[i*NDOF3+0]+SSA_values[i*NDOF2+0];
    10115                 vy[i]       = FS_values[i*NDOF3+1]+SSA_values[i*NDOF2+1];
    10116                 vzFS[i]     = FS_values[i*NDOF3+2];
    10117                 pressure[i] = FS_values[NUMVERTICES*NDOF3+i]*FSreconditioning;
    10118 
    10119                 /*Check solution*/
    10120                 if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
    10121                 if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
    10122                 if(xIsNan<IssmDouble>(vzFS[i]))     _error_("NaN found in solution vector");
    10123                 if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
    10124         }
    10125 
    10126         /*Get Vz*/
    10127         Input* vzSSA_input=inputs->GetInput(VzSSAEnum);
    10128         if (vzSSA_input){
    10129                 if (vzSSA_input->ObjectEnum()!=PentaInputEnum){
    10130                         _error_("Cannot compute Vel as VzSSA is of type " << EnumToStringx(vzSSA_input->ObjectEnum()));
    10131                 }
    10132                 GetInputListOnVertices(&vzSSA[0],VzSSAEnum);
    10133         }
    10134         else{
    10135                 _error_("Cannot update solution as VzSSA is not present");
    10136         }
    10137 
    10138         /*Now Compute vel*/
    10139         for(i=0;i<NUMVERTICES;i++) {
    10140                 vz[i]  = vzSSA[i]+vzFS[i];
    10141                 vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    10142         }
    10143 
    10144         /*Now, we have to move the previous Vx and Vy inputs  to old
    10145          * status, otherwise, we'll wipe them off: */
    10146         this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
    10147         this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
    10148         this->inputs->ChangeEnum(VzEnum,VzPicardEnum);
    10149         this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
    10150 
    10151         /*Add vx and vy as inputs to the tria element: */
    10152         this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
    10153         this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
    10154         this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
    10155         this->inputs->AddInput(new PentaInput(VzFSEnum,vzFS,P1Enum));
    10156         this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
    10157         this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    10158 
    10159         /*Free ressources:*/
    10160         xDelete<int>(doflistSSA);
    10161         xDelete<int>(doflistFSv);
    10162         xDelete<int>(doflistFSp);
    10163         xDelete<int>(cs_list);
    10164 }
    10165 /*}}}*/
    10166 /*FUNCTION Penta::InputUpdateFromSolutionStressbalanceL1L2 {{{*/
    10167 void  Penta::InputUpdateFromSolutionStressbalanceL1L2(IssmDouble* solution){
    10168 
    10169         const int    numdof=NDOF2*NUMVERTICES;
    10170 
    10171         int     i;
    10172         IssmDouble  rho_ice,g;
    10173         IssmDouble  values[numdof];
    10174         IssmDouble  vx[NUMVERTICES];
    10175         IssmDouble  vy[NUMVERTICES];
    10176         IssmDouble  vz[NUMVERTICES];
    10177         IssmDouble  vel[NUMVERTICES];
    10178         IssmDouble  pressure[NUMVERTICES];
    10179         IssmDouble  surface[NUMVERTICES];
    10180         IssmDouble  xyz_list[NUMVERTICES][3];
    10181         int    *doflist = NULL;
    10182         Penta  *penta   = NULL;
    10183 
    10184         /*Get dof list: */
    10185         GetDofList(&doflist,L1L2ApproximationEnum,GsetEnum);
    10186 
    10187         /*Use the dof list to index into the solution vector: */
    10188         for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    10189 
    10190         /*Transform solution in Cartesian Space*/
    10191         ::TransformSolutionCoord(&values[0],nodes,NUMVERTICES2D,XYEnum); /*2D: only the first 3 nodes are taken*/
    10192 
    10193         /*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */
    10194         for(i=0;i<3;i++){
    10195                 vx[i]  =values[i*NDOF2+0];
    10196                 vy[i]  =values[i*NDOF2+1];
    10197                 vx[i+3]=vx[i];
    10198                 vy[i+3]=vy[i];
    10199 
    10200                 /*Check solution*/
    10201                 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    10202                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
    10203         }
    10204 
    10205         /*Get parameters fro pressure computation*/
    10206         rho_ice=matpar->GetRhoIce();
    10207         g=matpar->GetG();
    10208 
    10209         /*Start looping over all elements above current element and update all inputs*/
    10210         penta=this;
    10211         for(;;){
    10212 
    10213                 /*Get node data: */
    10214                 ::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
    10215 
    10216                 /*Now Compute vel*/
    10217                 GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
    10218                 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);
    10219 
    10220                 /*Now compute pressure*/
    10221                 GetInputListOnVertices(&surface[0],SurfaceEnum);
    10222                 for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
    10223 
    10224                 /*Now, we have to move the previous Vx and Vy inputs  to old
    10225                  * status, otherwise, we'll wipe them off: */
    10226                 penta->inputs->ChangeEnum(VxEnum,VxPicardEnum);
    10227                 penta->inputs->ChangeEnum(VyEnum,VyPicardEnum);
    10228                 penta->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
    10229 
    10230                 /*Add vx and vy as inputs to the tria element: */
    10231                 penta->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
    10232                 penta->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
    10233                 penta->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
    10234                 penta->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    10235 
    10236                 /*Stop if we have reached the surface*/
    10237                 if (penta->IsOnSurface()) break;
    10238 
    10239                 /* get upper Penta*/
    10240                 penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
    10241         }
    10242 
    10243         /*Free ressources:*/
    10244         xDelete<int>(doflist);
    10245 }
    10246 /*}}}*/
    10247 /*FUNCTION Penta::InputUpdateFromSolutionStressbalanceHO {{{*/
    10248 void  Penta::InputUpdateFromSolutionStressbalanceHO(IssmDouble* solution){
    10249 
    10250         int         i;
    10251         IssmDouble  rho_ice,g;
    10252         IssmDouble  xyz_list[NUMVERTICES][3];
    10253         int        *doflist = NULL;
    10254 
    10255         /*Fetch number of nodes and dof for this finite element*/
    10256         int numnodes = this->NumberofNodes();
    10257         int numdof   = numnodes*NDOF2;
    10258 
    10259         /*Fetch dof list and allocate solution vectors*/
    10260         GetDofList(&doflist,HOApproximationEnum,GsetEnum);
    10261         IssmDouble* values    = xNew<IssmDouble>(numdof);
    10262         IssmDouble* vx        = xNew<IssmDouble>(numnodes);
    10263         IssmDouble* vy        = xNew<IssmDouble>(numnodes);
    10264         IssmDouble* vz        = xNew<IssmDouble>(numnodes);
    10265         IssmDouble* vel       = xNew<IssmDouble>(numnodes);
    10266         IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
    10267         IssmDouble* surface   = xNew<IssmDouble>(numnodes);
    10268 
    10269         /*Use the dof list to index into the solution vector: */
    10270         for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    10271 
    10272         /*Transform solution in Cartesian Space*/
    10273         ::TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum);
    10274         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    10275 
    10276         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    10277         for(i=0;i<numnodes;i++){
    10278                 vx[i]=values[i*NDOF2+0];
    10279                 vy[i]=values[i*NDOF2+1];
    10280 
    10281                 /*Check solution*/
    10282                 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    10283                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
    10284         }
    10285 
    10286         /*Get Vz and compute vel*/
    10287         GetInputListOnNodes(&vz[0],VzEnum,0.);
    10288         for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    10289 
    10290         /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D,
    10291          *so the pressure is just the pressure at the z elevation: */
    10292         rho_ice=matpar->GetRhoIce();
    10293         g=matpar->GetG();
    10294         GetInputListOnVertices(&surface[0],SurfaceEnum);
    10295         for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
    10296 
    10297         /*Now, we have to move the previous Vx and Vy inputs  to old
    10298          * status, otherwise, we'll wipe them off: */
    10299         this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
    10300         this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
    10301         this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
    10302 
    10303         /*Add vx and vy as inputs to the tria element: */
    10304         this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
    10305         this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
    10306         this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
    10307         this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    10308 
    10309         /*Free ressources:*/
    10310         xDelete<IssmDouble>(surface);
    10311         xDelete<IssmDouble>(pressure);
    10312         xDelete<IssmDouble>(vel);
    10313         xDelete<IssmDouble>(vz);
    10314         xDelete<IssmDouble>(vy);
    10315         xDelete<IssmDouble>(vx);
    10316         xDelete<IssmDouble>(values);
    10317         xDelete<int>(doflist);
    10318 }
    10319 /*}}}*/
    10320 /*FUNCTION Penta::InputUpdateFromSolutionStressbalanceHOFS {{{*/
    10321 void  Penta::InputUpdateFromSolutionStressbalanceHOFS(IssmDouble* solution){
    10322 
    10323         const int    numdofHO  = NDOF2*NUMVERTICES;
    10324         const int    numdofFSv = NDOF3*NUMVERTICES;
    10325         const int    numdofFSp = NDOF1*NUMVERTICES;
    10326 
    10327         int        i;
    10328         IssmDouble HO_values[numdofHO];
    10329         IssmDouble FS_values[numdofFSv+numdofFSp];
    10330         IssmDouble vx[NUMVERTICES];
    10331         IssmDouble vy[NUMVERTICES];
    10332         IssmDouble vz[NUMVERTICES];
    10333         IssmDouble vzHO[NUMVERTICES];
    10334         IssmDouble vzFS[NUMVERTICES];
    10335         IssmDouble vel[NUMVERTICES];
    10336         IssmDouble pressure[NUMVERTICES];
    10337         IssmDouble xyz_list[NUMVERTICES][3];
    10338         IssmDouble FSreconditioning;
    10339         int*       doflistHO        = NULL;
    10340         int*       doflistFSv        = NULL;
    10341         int*       doflistFSp = NULL;
    10342 
    10343         /*Prepare coordinate system list*/
    10344         int* cs_list = xNew<int>(2*NUMVERTICES);
    10345         for(i=0;i<NUMVERTICES;i++) cs_list[i]             = XYZEnum;
    10346         for(i=0;i<NUMVERTICES;i++) cs_list[NUMVERTICES+i] = PressureEnum;
    10347 
    10348         /*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
    10349         GetDofList(&doflistHO,HOApproximationEnum,GsetEnum);
    10350         GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
    10351         GetDofListPressure(&doflistFSp,GsetEnum);
    10352         this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    10353 
    10354         /*Get node data: */
    10355         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    10356 
    10357         /*Use the dof list to index into the solution vector: */
    10358         for(i=0;i<numdofHO;i++)  HO_values[i]=solution[doflistHO[i]];
    10359         for(i=0;i<numdofFSv;i++) FS_values[i]=solution[doflistFSv[i]];
    10360         for(i=0;i<numdofFSp;i++) FS_values[numdofFSv+i]=solution[doflistFSp[i]];
    10361 
    10362         /*Transform solution in Cartesian Space*/
    10363         ::TransformSolutionCoord(&HO_values[0],this->nodes,NUMVERTICES,XYEnum);
    10364         ::TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
    10365 
    10366         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    10367         for(i=0;i<NUMVERTICES;i++){
    10368                 vx[i]       = FS_values[i*NDOF3+0]+HO_values[i*NDOF2+0];
    10369                 vy[i]       = FS_values[i*NDOF3+1]+HO_values[i*NDOF2+1];
    10370                 vzFS[i]     = FS_values[i*NDOF3+2];
    10371                 pressure[i] = FS_values[NUMVERTICES*NDOF3+i]*FSreconditioning;
    10372 
    10373                 /*Check solution*/
    10374                 if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
    10375                 if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
    10376                 if(xIsNan<IssmDouble>(vzFS[i]))     _error_("NaN found in solution vector");
    10377                 if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
    10378         }
    10379 
    10380         /*Get Vz*/
    10381         Input* vzHO_input=inputs->GetInput(VzHOEnum);
    10382         if (vzHO_input){
    10383                 if (vzHO_input->ObjectEnum()!=PentaInputEnum){
    10384                         _error_("Cannot compute Vel as VzHO is of type " << EnumToStringx(vzHO_input->ObjectEnum()));
    10385                 }
    10386                 GetInputListOnVertices(&vzHO[0],VzHOEnum);
    10387         }
    10388         else{
    10389                 _error_("Cannot update solution as VzHO is not present");
    10390         }
    10391 
    10392         /*Now Compute vel*/
    10393         for(i=0;i<NUMVERTICES;i++) {
    10394                 vz[i]  = vzHO[i]+vzFS[i];
    10395                 vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    10396         }
    10397 
    10398         /*Now, we have to move the previous Vx and Vy inputs  to old
    10399          * status, otherwise, we'll wipe them off: */
    10400         this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
    10401         this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
    10402         this->inputs->ChangeEnum(VzEnum,VzPicardEnum);
    10403         this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
    10404 
    10405         /*Add vx and vy as inputs to the tria element: */
    10406         this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
    10407         this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
    10408         this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
    10409         this->inputs->AddInput(new PentaInput(VzFSEnum,vzFS,P1Enum));
    10410         this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
    10411         this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    10412 
    10413         /*Free ressources:*/
    10414         xDelete<int>(doflistHO);
    10415         xDelete<int>(doflistFSv);
    10416         xDelete<int>(doflistFSp);
    10417         xDelete<int>(cs_list);
    10418 }
    10419 /*}}}*/
    10420 /*FUNCTION Penta::InputUpdateFromSolutionStressbalanceSIA {{{*/
    10421 void  Penta::InputUpdateFromSolutionStressbalanceSIA(IssmDouble* solution){
    10422 
    10423         int         numnodes = this->NumberofNodes();
    10424         int         numdof=NDOF2*numnodes;
    10425 
    10426         int     i;
    10427         IssmDouble  rho_ice,g;
    10428         IssmDouble  values[numdof];
    10429         IssmDouble  vx[NUMVERTICES];
    10430         IssmDouble  vy[NUMVERTICES];
    10431         IssmDouble  vz[NUMVERTICES];
    10432         IssmDouble  vel[NUMVERTICES];
    10433         IssmDouble  pressure[NUMVERTICES];
    10434         IssmDouble  surface[NUMVERTICES];
    10435         IssmDouble  xyz_list[NUMVERTICES][3];
    10436         int*    doflist = NULL;
    10437 
    10438         /*Get dof list: */
    10439         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    10440 
    10441         /*Get node data: */
    10442         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    10443 
    10444         /*Use the dof list to index into the solution vector: */
    10445         for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    10446 
    10447         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    10448         for(i=0;i<NUMVERTICES;i++){
    10449                 vx[i]=values[i*NDOF2+0];
    10450                 vy[i]=values[i*NDOF2+1];
    10451 
    10452                 /*Check solution*/
    10453                 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    10454                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
    10455         }
    10456 
    10457         /*Now Compute vel*/
    10458         GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
    10459         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);
    10460 
    10461         /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D,
    10462          *so the pressure is just the pressure at the z elevation: */
    10463         rho_ice=matpar->GetRhoIce();
    10464         g=matpar->GetG();
    10465         GetInputListOnVertices(&surface[0],SurfaceEnum);
    10466         for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
    10467 
    10468         /*Now, we have to move the previous Vx and Vy inputs  to old
    10469          * status, otherwise, we'll wipe them off: */
    10470         this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
    10471         this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
    10472         this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
    10473 
    10474         /*Add vx and vy as inputs to the tria element: */
    10475         this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
    10476         this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
    10477         this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
    10478         this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    10479 
    10480         /*Free ressources:*/
    10481         xDelete<int>(doflist);
    10482 }
    10483 /*}}}*/
    10484 /*FUNCTION Penta::InputUpdateFromSolutionStressbalanceVert {{{*/
    10485 void  Penta::InputUpdateFromSolutionStressbalanceVert(IssmDouble* solution){
    10486 
    10487         int          numnodes = this->NumberofNodes();
    10488         int          numdof=NDOF1*numnodes;
    10489 
    10490         int          i;
    10491         int          approximation;
    10492         IssmDouble   rho_ice,g;
    10493         IssmDouble   values[numdof];
    10494         IssmDouble   vx[NUMVERTICES];
    10495         IssmDouble   vy[NUMVERTICES];
    10496         IssmDouble   vz[NUMVERTICES];
    10497         IssmDouble   vzSSA[NUMVERTICES];
    10498         IssmDouble   vzHO[NUMVERTICES];
    10499         IssmDouble   vzFS[NUMVERTICES];
    10500         IssmDouble   vel[NUMVERTICES];
    10501         IssmDouble   pressure[NUMVERTICES];
    10502         IssmDouble   surface[NUMVERTICES];
    10503         IssmDouble   xyz_list[NUMVERTICES][3];
    10504         int*         doflist      = NULL;
    10505 
    10506         /*Get the approximation and do nothing if the element in FS or None*/
    10507         inputs->GetInputValue(&approximation,ApproximationEnum);
    10508         if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
    10509                 return;
    10510         }
    10511 
    10512         /*Get dof list and vertices coordinates: */
    10513         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    10514         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    10515 
    10516         /*Use the dof list to index into the solution vector vz: */
    10517         for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    10518         for(i=0;i<NUMVERTICES;i++){
    10519                 vz[i]=values[i*NDOF1+0];
    10520 
    10521                 /*Check solution*/
    10522                 if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
    10523         }
    10524 
    10525         /*Get Vx and Vy*/
    10526         GetInputListOnVertices(&vx[0],VxEnum,0.0); //default is 0
    10527         GetInputListOnVertices(&vy[0],VyEnum,0.0); //default is 0
    10528 
    10529         /*Do some modifications if we actually have a HOFS or SSAFS element*/
    10530         if(approximation==HOFSApproximationEnum){
    10531                 Input* vzFS_input=inputs->GetInput(VzFSEnum);
    10532                 if (vzFS_input){
    10533                         if (vzFS_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
    10534                         GetInputListOnVertices(&vzFS[0],VzFSEnum);
    10535                 }
    10536                 else _error_("Cannot compute Vz as VzFS in not present in HOFS element");
    10537                 for(i=0;i<NUMVERTICES;i++){
    10538                         vzHO[i]=vz[i];
    10539                         vz[i]=vzHO[i]+vzFS[i];
    10540                 }
    10541         }
    10542         else if(approximation==SSAFSApproximationEnum){
    10543                 Input* vzFS_input=inputs->GetInput(VzFSEnum);
    10544                 if (vzFS_input){
    10545                         if (vzFS_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
    10546                         GetInputListOnVertices(&vzFS[0],VzFSEnum);
    10547                 }
    10548                 else _error_("Cannot compute Vz as VzFS in not present in SSAFS element");
    10549                 for(i=0;i<NUMVERTICES;i++){
    10550                         vzSSA[i]=vz[i];
    10551                         vz[i]=vzSSA[i]+vzFS[i];
    10552                 }
    10553         }
    10554 
    10555         /*Now Compute vel*/
    10556         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);
    10557 
    10558         /*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D,
    10559          *so the pressure is just the pressure at the z elevation: except it this is a HOFS element */
    10560         if(approximation!=HOFSApproximationEnum &&  approximation!=SSAFSApproximationEnum){
    10561                 rho_ice=matpar->GetRhoIce();
    10562                 g=matpar->GetG();
    10563                 GetInputListOnVertices(&surface[0],SurfaceEnum);
    10564                 for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
    10565         }
    10566 
    10567         /*Now, we have to move the previous Vz inputs to old
    10568          * status, otherwise, we'll wipe them off and add the new inputs: */
    10569         this->inputs->ChangeEnum(VzEnum,VzPicardEnum);
    10570 
    10571         if(approximation!=HOFSApproximationEnum && approximation!=SSAFSApproximationEnum){
    10572                 this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
    10573                 this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    10574         }
    10575         else if(approximation==HOFSApproximationEnum){
    10576                 this->inputs->AddInput(new PentaInput(VzHOEnum,vzHO,P1Enum));
    10577         }
    10578         else if(approximation==SSAFSApproximationEnum){
    10579                 this->inputs->AddInput(new PentaInput(VzSSAEnum,vzSSA,P1Enum));
    10580         }
    10581         this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
    10582         this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
    10583 
    10584         /*Free ressources:*/
    10585         xDelete<int>(doflist);
    10586 }
    10587 /*}}}*/
    10588 /*FUNCTION Penta::InputUpdateFromSolutionStressbalanceFS {{{*/
    10589 void  Penta::InputUpdateFromSolutionStressbalanceFS(IssmDouble* solution){
    10590 
    10591         int          i;
    10592         int*         vdoflist=NULL;
    10593         int*         pdoflist=NULL;
    10594         IssmDouble   FSreconditioning;
    10595 
    10596         /*Fetch number of nodes and dof for this finite element*/
    10597         int vnumnodes = this->NumberofNodesVelocity();
    10598         int pnumnodes = this->NumberofNodesPressure();
    10599         int vnumdof   = vnumnodes*NDOF3;
    10600         int pnumdof   = pnumnodes*NDOF1;
    10601 
    10602         /*Initialize values*/
    10603         IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
    10604         IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
    10605         IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
    10606         IssmDouble* vz       = xNew<IssmDouble>(vnumnodes);
    10607         IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
    10608         IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
    10609 
    10610         /*Prepare coordinate system list*/
    10611         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    10612         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
    10613         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    10614 
    10615         /*Get dof list: */
    10616         GetDofListVelocity(&vdoflist,GsetEnum);
    10617         GetDofListPressure(&pdoflist,GsetEnum);
    10618 
    10619         /*Use the dof list to index into the solution vector: */
    10620         for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
    10621         for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
    10622 
    10623         /*Transform solution in Cartesian Space*/
    10624         ::TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
    10625 
    10626         /*Ok, we have vx and vy in values, fill in all arrays: */
    10627         for(i=0;i<vnumnodes;i++){
    10628                 vx[i] = values[i*NDOF3+0];
    10629                 vy[i] = values[i*NDOF3+1];
    10630                 vz[i] = values[i*NDOF3+2];
    10631                 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    10632                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
    10633                 if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
    10634         }
    10635         for(i=0;i<pnumnodes;i++){
    10636                 pressure[i] = values[vnumdof+i];
    10637                 if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
    10638         }
    10639 
    10640         /*Recondition pressure and compute vel: */
    10641         this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    10642         for(i = 0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
    10643         for(i = 0;i<vnumnodes;i++) vel[i]      = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    10644 
    10645         /*Now, we have to move the previous inputs  to old
    10646          * status, otherwise, we'll wipe them off: */
    10647         this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
    10648         this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
    10649         this->inputs->ChangeEnum(VzEnum,VzPicardEnum);
    10650         this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
    10651 
    10652         /*Add vx and vy as inputs to the tria element: */
    10653         this->inputs->AddInput(new PentaInput(VxEnum,vx,P1Enum));
    10654         this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
    10655         this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
    10656         this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
    10657         this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
    10658 
    10659         /*Free ressources:*/
    10660         xDelete<IssmDouble>(pressure);
    10661         xDelete<IssmDouble>(vel);
    10662         xDelete<IssmDouble>(vz);
    10663         xDelete<IssmDouble>(vy);
    10664         xDelete<IssmDouble>(vx);
    10665         xDelete<IssmDouble>(values);
    10666         xDelete<int>(vdoflist);
    10667         xDelete<int>(pdoflist);
    10668         xDelete<int>(cs_list);
    106699331}
    106709332/*}}}*/
    … …  
    109529614}
    109539615/*}}}*/
    10954 /*FUNCTION Penta::InputUpdateFromSolutionHydrologyDCInefficient{{{*/
    10955 void  Penta::InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution){
    10956 
    10957         const int   numdof   = NDOF1*NUMVERTICES;
    10958         const int   numdof2d = NDOF1*NUMVERTICES2D;
    10959         int*        doflist  = NULL;
    10960         bool        converged;
    10961         IssmDouble  values[numdof];
    10962         IssmDouble  residual[numdof];
    10963         IssmDouble  penalty_factor;
    10964         IssmDouble  kmax, kappa, h_max;
    10965         Penta      *penta    = NULL;
    10966 
    10967         /*If not on bed, return*/
    10968         if (!IsOnBed()) return;
    10969 
    10970         /*Get dof list: */
    10971         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    10972 
    10973         /*Use the dof list to index into the solution vector and extrude it */
    10974         for(int i=0;i<numdof2d;i++){
    10975                 values[i]         =solution[doflist[i]];
    10976                 values[i+numdof2d]=values[i];
    10977                 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
    10978         }
    10979 
    10980         /*If converged keep the residual in mind*/
    10981         this->inputs->GetInputValue(&converged,ConvergedEnum);
    10982 
    10983         /*Get inputs*/
    10984         if(converged){
    10985                 this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum);
    10986                 this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
    10987 
    10988                 kappa=kmax*pow(10.,penalty_factor);
    10989 
    10990                 Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.     
    10991                 for(int i=0;i<NUMVERTICES2D;i++){
    10992                         tria->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
    10993                         if(values[i]>h_max){
    10994                                 residual[i]=kappa*(values[i]-h_max);
    10995                                 residual[i+numdof2d]=residual[i];
    10996                         }
    10997                         else{
    10998                                 residual[i]=0.0;
    10999                                 residual[i+numdof2d]=residual[i];
    11000                         }
    11001                 }
    11002                 delete tria->material; delete tria;
    11003         }
    11004 
    11005         /*Start looping over all elements above current element and update all inputs*/
    11006         penta=this;
    11007         for(;;){
    11008                 /*Add input to the element: */
    11009                 penta->inputs->AddInput(new PentaInput(SedimentHeadEnum,values,P1Enum));
    11010                 penta->inputs->AddInput(new PentaInput(SedimentHeadResidualEnum,residual,P1Enum));
    11011 
    11012                 /*Stop if we have reached the surface*/
    11013                 if (penta->IsOnSurface()) break;
    11014 
    11015                 /* get upper Penta*/
    11016                 penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
    11017         }
    11018 
    11019         /*Free ressources:*/
    11020         xDelete<int>(doflist);
    11021 }
    11022 /*}}}*/
    110239616/*FUNCTION Penta::UpdateConstraintsL2ProjectionEPL{{{*/
    110249617void  Penta::UpdateConstraintsL2ProjectionEPL(void){
Note: See TracChangeset for help on using the changeset viewer.