Changeset 16783


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

NEW: InputUpdateFromSolution is now entirely done by analyses

Location:
issm/trunk-jpl/src/c
Files:
16 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){

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

    r16778 r16783  
    6060                void  InputUpdateFromConstant(IssmDouble constant, int name);
    6161                void  InputUpdateFromConstant(int constant, int name);
    62                 void  InputUpdateFromSolution(IssmDouble* solutiong);
    6362                void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
    6463                #ifdef _HAVE_DAKOTA_
    … …  
    232231                void           InputChangeName(int input_enum, int enum_type_old);
    233232                void             InputExtrude(int enum_type,int object_type);
    234                 void           InputUpdateFromSolutionMasstransport(IssmDouble* solutiong);
    235                 void           InputUpdateFromSolutionFreeSurfaceTop(IssmDouble* solutiong);
    236                 void           InputUpdateFromSolutionFreeSurfaceBase(IssmDouble* solutiong);
    237233                void           InputUpdateFromSolutionOneDof(IssmDouble* solutiong,int enum_type);
    238234                void           InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solutiong,int enum_type);
    … …  
    290286                ElementMatrix* CreateJacobianStressbalanceHO(void);
    291287                ElementMatrix* CreateJacobianStressbalanceFS(void);
    292                 void           InputUpdateFromSolutionStressbalanceHoriz( IssmDouble* solutiong);
    293                 void           InputUpdateFromSolutionStressbalanceSSA( IssmDouble* solutiong);
    294                 void           InputUpdateFromSolutionStressbalanceSSAHO( IssmDouble* solutiong);
    295                 void           InputUpdateFromSolutionStressbalanceSSAFS( IssmDouble* solutiong);
    296                 void           InputUpdateFromSolutionStressbalanceL1L2( IssmDouble* solutiong);
    297                 void           InputUpdateFromSolutionStressbalanceHO( IssmDouble* solutiong);
    298                 void           InputUpdateFromSolutionStressbalanceHOFS( IssmDouble* solutiong);
    299                 void           InputUpdateFromSolutionStressbalanceSIA( IssmDouble* solutiong);
    300                 void           InputUpdateFromSolutionStressbalanceVert( IssmDouble* solutiong);
    301                 void           InputUpdateFromSolutionStressbalanceFS( IssmDouble* solutiong);
    302288                ElementVector* CreatePVectorCouplingSSAFS(void);
    303289                ElementVector* CreatePVectorCouplingSSAFSViscous(void);
    … …  
    335321                ElementVector* CreatePVectorAdjointHO(void);
    336322                ElementVector* CreatePVectorAdjointFS(void);
    337                 void           InputUpdateFromSolutionAdjointHoriz( IssmDouble* solutiong);
    338                 void           InputUpdateFromSolutionAdjointFS( IssmDouble* solutiong);
    339323                #endif
    340324
    … …  
    351335                void           HydrologyEPLGetActive(Vector<IssmDouble>* active_vec);
    352336                void           HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
    353                 void           InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
    354337                void           ComputeEPLThickness(void);
    355338                void           UpdateConstraintsL2ProjectionEPL(void);
    … …  
    375358                ElementVector* CreatePVectorThermalShelf(void);
    376359                ElementVector* CreatePVectorThermalSheet(void);
    377                 void           InputUpdateFromSolutionThermal(IssmDouble* solutiong);
    378                 void           InputUpdateFromSolutionEnthalpy(IssmDouble* solutiong);
    379360                void           UpdateBasalConstraintsEnthalpy(void);
    380361                void           ComputeBasalMeltingrate(void);

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r16778 r16783  
    5454                /*}}}*/
    5555                /*Update virtual functions resolution: {{{*/
    56                 void  InputUpdateFromSolution(IssmDouble* solution){_error_("not implemented yet");};
    5756                void  InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int inputenum){_error_("not implemented yet");};
    5857                void  InputUpdateFromSolutionOneDof(IssmDouble* solution,int inputenum){_error_("not implemented yet");};

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

    r16778 r16783  
    18091809}
    18101810/*}}}*/
    1811 /*FUNCTION Tria::InputUpdateFromSolution {{{*/
    1812 void  Tria::InputUpdateFromSolution(IssmDouble* solution){
    1813 
    1814         /*retrive parameters: */
    1815         int analysis_type,extrusioninput;
    1816         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    1817 
    1818         /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
    1819         switch(analysis_type){
    1820                 #ifdef _HAVE_STRESSBALANCE_
    1821                 case StressbalanceAnalysisEnum:
    1822                         int approximation;
    1823                         inputs->GetInputValue(&approximation,ApproximationEnum);
    1824                         if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
    1825                                 InputUpdateFromSolutionStressbalanceFS(solution);
    1826                         }
    1827                         else if (approximation==SSAApproximationEnum || approximation==SIAApproximationEnum){
    1828                                 InputUpdateFromSolutionStressbalanceHoriz(solution);
    1829                         }
    1830                         else{
    1831                                 _error_("approximation not supported yet");
    1832                         }
    1833                         break;
    1834                 case StressbalanceSIAAnalysisEnum:
    1835                         InputUpdateFromSolutionStressbalanceHoriz(solution);
    1836                         break;
    1837                 #endif
    1838                 #ifdef _HAVE_CONTROL_
    1839                 case AdjointHorizAnalysisEnum:
    1840                         InputUpdateFromSolutionAdjointHoriz(solution);
    1841                         break;
    1842                 case AdjointBalancethicknessAnalysisEnum:
    1843                         InputUpdateFromSolutionOneDof(solution,AdjointEnum);
    1844                         break;
    1845                 #endif
    1846                 #ifdef _HAVE_HYDROLOGY_
    1847                 case HydrologyShreveAnalysisEnum:
    1848                         InputUpdateFromSolutionHydrologyShreve(solution);
    1849                         break;
    1850                 case HydrologyDCInefficientAnalysisEnum:
    1851                         InputUpdateFromSolutionHydrologyDCInefficient(solution);
    1852                         break;
    1853                 case HydrologyDCEfficientAnalysisEnum:
    1854                         InputUpdateFromSolutionOneDof(solution,EplHeadEnum);
    1855                         break;
    1856                 case L2ProjectionEPLAnalysisEnum:
    1857                         this->parameters->FindParam(&extrusioninput,InputToL2ProjectEnum);
    1858                         InputUpdateFromSolutionOneDof(solution,extrusioninput);
    1859                         break;
    1860                 #endif
    1861                 #ifdef _HAVE_DAMAGE_
    1862                 case DamageEvolutionAnalysisEnum:
    1863                         InputUpdateFromSolutionDamageEvolution(solution);
    1864                         break;
    1865                 #endif
    1866                 #ifdef _HAVE_BALANCED_
    1867                 case BalancethicknessAnalysisEnum:
    1868                         InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
    1869                         break;
    1870                 case BalancevelocityAnalysisEnum:
    1871                         InputUpdateFromSolutionOneDof(solution,VelEnum);
    1872                         break;
    1873                 case SmoothedSurfaceSlopeXAnalysisEnum:
    1874                         InputUpdateFromSolutionOneDof(solution,SurfaceSlopeXEnum);
    1875                         break;
    1876                 case SmoothedSurfaceSlopeYAnalysisEnum:
    1877                         InputUpdateFromSolutionOneDof(solution,SurfaceSlopeYEnum);
    1878                         break;
    1879                 #endif
    1880                 case L2ProjectionBaseAnalysisEnum:
    1881                         this->parameters->FindParam(&extrusioninput,InputToL2ProjectEnum);
    1882                         InputUpdateFromSolutionOneDof(solution,extrusioninput);
    1883                         break;
    1884                 case MasstransportAnalysisEnum:
    1885                         InputUpdateFromSolutionMasstransport(solution);
    1886                         break;
    1887                 case FreeSurfaceTopAnalysisEnum:
    1888                         InputUpdateFromSolutionOneDof(solution,SurfaceEnum);
    1889                         break;
    1890                 case FreeSurfaceBaseAnalysisEnum:
    1891                         InputUpdateFromSolutionOneDof(solution,BedEnum);
    1892                         break;
    1893                 case ExtrudeFromBaseAnalysisEnum: case ExtrudeFromTopAnalysisEnum:
    1894                         this->parameters->FindParam(&extrusioninput,InputToExtrudeEnum);
    1895                         InputUpdateFromSolutionOneDof(solution,extrusioninput);
    1896                         break;
    1897                 default:
    1898                         _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    1899         }
    1900 }
    1901 /*}}}*/
    19021811/*FUNCTION Tria::InputUpdateFromSolutionOneDof{{{*/
    19031812void  Tria::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){
    … …  
    19241833        /*Free ressources:*/
    19251834        xDelete<IssmDouble>(values);
    1926         xDelete<int>(doflist);
    1927 }
    1928 /*}}}*/
    1929 /*FUNCTION Tria::InputUpdateFromSolutionMasstransport{{{*/
    1930 void  Tria::InputUpdateFromSolutionMasstransport(IssmDouble* solution){
    1931 
    1932         /*Intermediaries*/
    1933         int        i,hydroadjustment;
    1934         int*       doflist=NULL;
    1935         IssmDouble rho_ice,rho_water,minthickness;
    1936 
    1937         /*Fetch number of nodes for this finite element*/
    1938         int numnodes = this->NumberofNodes();
    1939 
    1940         /*Fetch dof list and allocate solution vector*/
    1941         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    1942         IssmDouble* newthickness = xNew<IssmDouble>(numnodes);
    1943         IssmDouble* newbed       = xNew<IssmDouble>(numnodes);
    1944         IssmDouble* newsurface   = xNew<IssmDouble>(numnodes);
    1945         IssmDouble* oldthickness = xNew<IssmDouble>(numnodes);
    1946         IssmDouble* oldbed       = xNew<IssmDouble>(numnodes);
    1947         IssmDouble* oldsurface   = xNew<IssmDouble>(numnodes);
    1948         IssmDouble* phi          = xNew<IssmDouble>(numnodes);
    1949 
    1950         /*Use the dof list to index into the solution vector: */
    1951         this->parameters->FindParam(&minthickness,MasstransportMinThicknessEnum);
    1952         for(i=0;i<numnodes;i++){
    1953                 newthickness[i]=solution[doflist[i]];
    1954                 if(xIsNan<IssmDouble>(newthickness[i])) _error_("NaN found in solution vector");
    1955                 /*Constrain thickness to be at least 1m*/
    1956                 if(newthickness[i]<minthickness) newthickness[i]=minthickness;
    1957         }
    1958 
    1959         /*Get previous bed, thickness and surface*/
    1960         GetInputListOnNodes(&oldbed[0],BedEnum);
    1961         GetInputListOnNodes(&oldsurface[0],SurfaceEnum);
    1962         GetInputListOnNodes(&oldthickness[0],ThicknessEnum);
    1963         GetInputListOnNodes(&phi[0],MaskGroundediceLevelsetEnum);
    1964 
    1965         /*Find MasstransportHydrostaticAdjustment to figure out how to update the geometry:*/
    1966         this->parameters->FindParam(&hydroadjustment,MasstransportHydrostaticAdjustmentEnum);
    1967         rho_ice=matpar->GetRhoIce();
    1968         rho_water=matpar->GetRhoWater();
    1969 
    1970         for(i=0;i<numnodes;i++) {
    1971                 /*If shelf: hydrostatic equilibrium*/
    1972                 if (phi[i]>0.){
    1973                         newsurface[i] = oldbed[i]+newthickness[i]; //surface = oldbed + newthickness
    1974                         newbed[i]     = oldbed[i];                 //same bed: do nothing
    1975                 }
    1976                 else{ //this is an ice shelf
    1977                         if(hydroadjustment==AbsoluteEnum){
    1978                                 newsurface[i] = newthickness[i]*(1-rho_ice/rho_water);
    1979                                 newbed[i]     = newthickness[i]*(-rho_ice/rho_water);
    1980                         }
    1981                         else if(hydroadjustment==IncrementalEnum){
    1982                                 newsurface[i] = oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH
    1983                                 newbed[i]     = oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed               = oldbed + di * dH
    1984                         }
    1985                         else _error_("Hydrostatic adjustment " << hydroadjustment << " (" << EnumToStringx(hydroadjustment) << ") not supported yet");
    1986                 }
    1987         }
    1988 
    1989         /*Add input to the element: */
    1990         this->inputs->AddInput(new TriaInput(ThicknessEnum,newthickness,P1Enum));
    1991         this->inputs->AddInput(new TriaInput(SurfaceEnum,newsurface,P1Enum));
    1992         this->inputs->AddInput(new TriaInput(BedEnum,newbed,P1Enum));
    1993 
    1994         /*Free ressources:*/
    1995         xDelete<IssmDouble>(newthickness);
    1996         xDelete<IssmDouble>(newbed);
    1997         xDelete<IssmDouble>(newsurface);
    1998         xDelete<IssmDouble>(oldthickness);
    1999         xDelete<IssmDouble>(oldbed);
    2000         xDelete<IssmDouble>(oldsurface);
    2001         xDelete<IssmDouble>(phi);
    20021835        xDelete<int>(doflist);
    20031836}
    … …  
    44144247
    44154248        return y;
    4416 }
    4417 /*}}}*/
    4418 /*FUNCTION Tria::InputUpdateFromSolutionStressbalanceHoriz {{{*/
    4419 void  Tria::InputUpdateFromSolutionStressbalanceHoriz(IssmDouble* solution){
    4420 
    4421         int        i;
    4422         IssmDouble rho_ice,g;
    4423         int*       doflist=NULL;
    4424 
    4425         /*Fetch number of nodes and dof for this finite element*/
    4426         int numnodes = this->NumberofNodes();
    4427         int numdof   = numnodes*NDOF2;
    4428 
    4429         /*Fetch dof list and allocate solution vectors*/
    4430         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    4431         IssmDouble* values    = xNew<IssmDouble>(numdof);
    4432         IssmDouble* vx        = xNew<IssmDouble>(numnodes);
    4433         IssmDouble* vy        = xNew<IssmDouble>(numnodes);
    4434         IssmDouble* vz        = xNew<IssmDouble>(numnodes);
    4435         IssmDouble* vel       = xNew<IssmDouble>(numnodes);
    4436         IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
    4437         IssmDouble* thickness = xNew<IssmDouble>(numnodes);
    4438 
    4439         /*Use the dof list to index into the solution vector: */
    4440         for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    4441 
    4442         /*Transform solution in Cartesian Space*/
    4443         ::TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
    4444 
    4445         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    4446         for(i=0;i<numnodes;i++){
    4447                 vx[i]=values[i*NDOF2+0];
    4448                 vy[i]=values[i*NDOF2+1];
    4449 
    4450                 /*Check solution*/
    4451                 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    4452                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
    4453         }
    4454 
    4455         /*Get Vz and compute vel*/
    4456         GetInputListOnNodes(&vz[0],VzEnum,0.);
    4457         for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    4458 
    4459         /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
    4460          *so the pressure is just the pressure at the bedrock: */
    4461         rho_ice=matpar->GetRhoIce();
    4462         g=matpar->GetG();
    4463         GetInputListOnNodes(&thickness[0],ThicknessEnum);
    4464         for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*thickness[i];
    4465 
    4466         /*Now, we have to move the previous Vx and Vy inputs  to old
    4467          * status, otherwise, we'll wipe them off: */
    4468         this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
    4469         this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
    4470         this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
    4471 
    4472         /*Add vx and vy as inputs to the tria element: */
    4473         this->inputs->AddInput(new TriaInput(VxEnum,vx,P1Enum));
    4474         this->inputs->AddInput(new TriaInput(VyEnum,vy,P1Enum));
    4475         this->inputs->AddInput(new TriaInput(VelEnum,vel,P1Enum));
    4476         this->inputs->AddInput(new TriaInput(PressureEnum,pressure,P1Enum));
    4477 
    4478         /*Free ressources:*/
    4479         xDelete<IssmDouble>(thickness);
    4480         xDelete<IssmDouble>(pressure);
    4481         xDelete<IssmDouble>(vel);
    4482         xDelete<IssmDouble>(vz);
    4483         xDelete<IssmDouble>(vy);
    4484         xDelete<IssmDouble>(vx);
    4485         xDelete<IssmDouble>(values);
    4486         xDelete<int>(doflist);
    4487 
    4488 }
    4489 /*}}}*/
    4490 /*FUNCTION Tria::InputUpdateFromSolutionStressbalanceFS {{{*/
    4491 void  Tria::InputUpdateFromSolutionStressbalanceFS(IssmDouble* solution){
    4492 
    4493         int          i;
    4494         int*         vdoflist=NULL;
    4495         int*         pdoflist=NULL;
    4496         IssmDouble   FSreconditioning;
    4497 
    4498         /*Fetch number of nodes and dof for this finite element*/
    4499         int vnumnodes = this->NumberofNodesVelocity();
    4500         int pnumnodes = this->NumberofNodesPressure();
    4501         int vnumdof   = vnumnodes*NDOF2;
    4502         int pnumdof   = pnumnodes*NDOF1;
    4503 
    4504         /*Initialize values*/
    4505         IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
    4506         IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
    4507         IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
    4508         IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
    4509         IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
    4510 
    4511         /*Prepare coordinate system list*/
    4512         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    4513         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
    4514         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    4515 
    4516         /*Get dof list: */
    4517         GetDofListVelocity(&vdoflist,GsetEnum);
    4518         GetDofListPressure(&pdoflist,GsetEnum);
    4519 
    4520         /*Use the dof list to index into the solution vector: */
    4521         for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
    4522         for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
    4523 
    4524         /*Transform solution in Cartesian Space*/
    4525         ::TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
    4526 
    4527         /*Ok, we have vx and vy in values, fill in all arrays: */
    4528         for(i=0;i<vnumnodes;i++){
    4529                 vx[i] = values[i*NDOF2+0];
    4530                 vy[i] = values[i*NDOF2+1];
    4531                 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    4532                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
    4533         }
    4534         for(i=0;i<pnumnodes;i++){
    4535                 pressure[i] = values[vnumdof+i];
    4536                 if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
    4537         }
    4538 
    4539         /*Recondition pressure and compute vel: */
    4540         this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    4541         for(i=0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
    4542         for(i=0;i<vnumnodes;i++) vel[i]      = sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
    4543 
    4544         /*Now, we have to move the previous inputs  to old
    4545          * status, otherwise, we'll wipe them off: */
    4546         this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
    4547         this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
    4548         this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
    4549 
    4550         /*Add vx and vy as inputs to the tria element: */
    4551         this->inputs->AddInput(new TriaInput(VxEnum,vx,P1Enum));
    4552         this->inputs->AddInput(new TriaInput(VyEnum,vy,P1Enum));
    4553         this->inputs->AddInput(new TriaInput(VelEnum,vel,P1Enum));
    4554         this->inputs->AddInput(new TriaInput(PressureEnum,pressure,P1Enum));
    4555 
    4556         /*Free ressources:*/
    4557         xDelete<IssmDouble>(pressure);
    4558         xDelete<IssmDouble>(vel);
    4559         xDelete<IssmDouble>(vy);
    4560         xDelete<IssmDouble>(vx);
    4561         xDelete<IssmDouble>(values);
    4562         xDelete<int>(vdoflist);
    4563         xDelete<int>(pdoflist);
    4564         xDelete<int>(cs_list);
    4565 }
    4566 /*}}}*/
    4567 /*FUNCTION Tria::InputUpdateFromSolutionStressbalanceSIA {{{*/
    4568 void  Tria::InputUpdateFromSolutionStressbalanceSIA(IssmDouble* solution){
    4569 
    4570         int        i;
    4571         IssmDouble rho_ice,g;
    4572         int*       doflist=NULL;
    4573 
    4574         /*Fetch number of nodes and dof for this finite element*/
    4575         int numnodes = this->NumberofNodes();
    4576         int numdof   = numnodes*NDOF2;
    4577 
    4578         /*Fetch dof list and allocate solution vectors*/
    4579         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    4580         IssmDouble* values    = xNew<IssmDouble>(numdof);
    4581         IssmDouble* vx        = xNew<IssmDouble>(numdof);
    4582         IssmDouble* vy        = xNew<IssmDouble>(numdof);
    4583         IssmDouble* vz        = xNew<IssmDouble>(numdof);
    4584         IssmDouble* vel       = xNew<IssmDouble>(numdof);
    4585         IssmDouble* pressure  = xNew<IssmDouble>(numdof);
    4586         IssmDouble* thickness = xNew<IssmDouble>(numdof);
    4587 
    4588         /*Use the dof list to index into the solution vector: */
    4589         for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    4590 
    4591         /*Transform solution in Cartesian Space*/
    4592         ::TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
    4593 
    4594         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    4595         for(i=0;i<numnodes;i++){
    4596                 vx[i]=values[i*NDOF2+0];
    4597                 vy[i]=values[i*NDOF2+1];
    4598 
    4599                 /*Check solution*/
    4600                 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    4601                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
    4602         }
    4603 
    4604         /*Get Vz and compute vel*/
    4605         GetInputListOnNodes(&vz[0],VzEnum,0.);
    4606         for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    4607 
    4608         /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
    4609          *so the pressure is just the pressure at the bedrock: */
    4610         rho_ice=matpar->GetRhoIce();
    4611         g=matpar->GetG();
    4612         GetInputListOnNodes(&thickness[0],ThicknessEnum);
    4613         for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*thickness[i];
    4614 
    4615         /*Now, we have to move the previous Vx and Vy inputs  to old
    4616          * status, otherwise, we'll wipe them off: */
    4617         this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
    4618         this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
    4619         this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
    4620 
    4621         /*Add vx and vy as inputs to the tria element: */
    4622         this->inputs->AddInput(new TriaInput(VxEnum,vx,P1Enum));
    4623         this->inputs->AddInput(new TriaInput(VyEnum,vy,P1Enum));
    4624         this->inputs->AddInput(new TriaInput(VelEnum,vel,P1Enum));
    4625         this->inputs->AddInput(new TriaInput(PressureEnum,pressure,P1Enum));
    4626 
    4627         /*Free ressources:*/
    4628         xDelete<IssmDouble>(thickness);
    4629         xDelete<IssmDouble>(pressure);
    4630         xDelete<IssmDouble>(vel);
    4631         xDelete<IssmDouble>(vz);
    4632         xDelete<IssmDouble>(vy);
    4633         xDelete<IssmDouble>(vx);
    4634         xDelete<IssmDouble>(values);
    4635         xDelete<int>(doflist);
    46364249}
    46374250/*}}}*/
    … …  
    63405953}
    63415954/*}}}*/
    6342 /*FUNCTION Tria::InputUpdateFromSolutionAdjointHoriz {{{*/
    6343 void  Tria::InputUpdateFromSolutionAdjointHoriz(IssmDouble* solution){
    6344 
    6345         int  i;
    6346         int* doflist=NULL;
    6347 
    6348         /*Fetch number of nodes and dof for this finite element*/
    6349         int numnodes = this->NumberofNodes();
    6350         int numdof   = numnodes*NDOF2;
    6351 
    6352         /*Fetch dof list and allocate solution vectors*/
    6353         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    6354         IssmDouble* values  = xNew<IssmDouble>(numdof);
    6355         IssmDouble* lambdax = xNew<IssmDouble>(numnodes);
    6356         IssmDouble* lambday = xNew<IssmDouble>(numnodes);
    6357 
    6358         /*Use the dof list to index into the solution vector: */
    6359         for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
    6360 
    6361         /*Transform solution in Cartesian Space*/
    6362         ::TransformSolutionCoord(&values[0],nodes,numnodes,XYEnum);
    6363 
    6364         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    6365         for(i=0;i<numnodes;i++){
    6366                 lambdax[i]=values[i*NDOF2+0];
    6367                 lambday[i]=values[i*NDOF2+1];
    6368 
    6369                 /*Check solution*/
    6370                 if(xIsNan<IssmDouble>(lambdax[i])) _error_("NaN found in solution vector");
    6371                 if(xIsNan<IssmDouble>(lambday[i])) _error_("NaN found in solution vector");
    6372         }
    6373 
    6374         /*Add vx and vy as inputs to the tria element: */
    6375         this->inputs->AddInput(new TriaInput(AdjointxEnum,lambdax,P1Enum));
    6376         this->inputs->AddInput(new TriaInput(AdjointyEnum,lambday,P1Enum));
    6377 
    6378         /*Free ressources:*/
    6379         xDelete<IssmDouble>(values);
    6380         xDelete<IssmDouble>(lambdax);
    6381         xDelete<IssmDouble>(lambday);
    6382         xDelete<int>(doflist);
    6383 }
    6384 /*}}}*/
    63855955/*FUNCTION Tria::GetVectorFromControlInputs{{{*/
    63865956void  Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){
    … …  
    71136683        xDelete<IssmDouble>(values);
    71146684        delete gauss;
    7115 }
    7116 /*}}}*/
    7117 /*FUNCTION Tria::InputUpdateFromSolutionHydrologyShreve{{{*/
    7118 void  Tria::InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution){
    7119 
    7120         /*Intermediary*/
    7121         int* doflist = NULL;
    7122 
    7123         /*Fetch number of nodes for this finite element*/
    7124         int numnodes = this->NumberofNodes();
    7125 
    7126         /*Fetch dof list and allocate solution vector*/
    7127         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    7128         IssmDouble* values    = xNew<IssmDouble>(numnodes);
    7129 
    7130         /*Use the dof list to index into the solution vector: */
    7131         for(int i=0;i<numnodes;i++){
    7132                 values[i]=solution[doflist[i]];
    7133                 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
    7134                 if (values[i]<10e-10) values[i]=10e-10; //correcting the water column to positive values
    7135         }
    7136 
    7137         /*Add input to the element: */
    7138         this->inputs->AddInput(new TriaInput(WatercolumnEnum,values,P1Enum));
    7139 
    7140         /*Free ressources:*/
    7141         xDelete<IssmDouble>(values);
    7142         xDelete<int>(doflist);
    7143 }
    7144 /*}}}*/
    7145 /*FUNCTION Tria::InputUpdateFromSolutionHydrologyDCInefficient{{{*/
    7146 void  Tria::InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution){
    7147 
    7148         /*Intermediaries*/
    7149         int        *doflist  = NULL;
    7150         bool        converged;
    7151         IssmDouble  penalty_factor, dt;
    7152         IssmDouble  kmax, kappa, h_max;
    7153 
    7154         /*Fetch number of nodes for this finite element*/
    7155         int numnodes = this->NumberofNodes();
    7156 
    7157         /*Fetch dof list and allocate solution vector*/
    7158         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    7159         IssmDouble* values   = xNew<IssmDouble>(numnodes);
    7160         IssmDouble* residual = xNew<IssmDouble>(numnodes);
    7161 
    7162         /*Use the dof list to index into the solution vector: */
    7163         for(int i=0;i<numnodes;i++){
    7164                 values[i]=solution[doflist[i]];
    7165                 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
    7166         }
    7167 
    7168         /*If converged keep the residual in mind*/
    7169         this->inputs->GetInputValue(&converged,ConvergedEnum);
    7170 
    7171         /*Get inputs*/
    7172         if(converged){
    7173                 this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum);
    7174                 this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
    7175 
    7176                 kappa=kmax*pow(10.,penalty_factor);
    7177 
    7178                 for(int i=0;i<numnodes;i++){
    7179                         this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
    7180                         if(values[i]>h_max){
    7181                                 residual[i]=kappa*(values[i]-h_max);
    7182                         }
    7183                         else{
    7184                                 residual[i]=0.0;
    7185                         }
    7186                 }
    7187                 this->inputs->AddInput(new TriaInput(SedimentHeadResidualEnum,residual,P1Enum));
    7188         }
    7189 
    7190         /*Add input to the element: */
    7191         this->inputs->AddInput(new TriaInput(SedimentHeadEnum,values,P1Enum));
    7192 
    7193         /*Free ressources:*/
    7194         xDelete<IssmDouble>(values);
    7195         xDelete<IssmDouble>(residual);
    7196         xDelete<int>(doflist);
    71976685}
    71986686/*}}}*/
    … …  
    85678055}
    85688056/*}}}*/
    8569 /*FUNCTION Tria::InputUpdateFromSolutionDamageEvolution {{{*/
    8570 void  Tria::InputUpdateFromSolutionDamageEvolution(IssmDouble* solution){
    8571 
    8572         const int    numdof=NDOF1*NUMVERTICES;
    8573 
    8574         int         i;
    8575         IssmDouble  values[numdof];
    8576         IssmDouble  max_damage;
    8577         int                     *doflist = NULL;
    8578 
    8579         /*Get dof list: */
    8580         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    8581 
    8582         /*Get user-supplied max_damage: */
    8583         this->parameters->FindParam(&max_damage,DamageMaxDamageEnum);
    8584 
    8585         /*Use the dof list to index into the solution vector: */
    8586         for(i=0;i<numdof;i++){
    8587                 values[i]=solution[doflist[i]];
    8588                 /*Check solution*/
    8589                 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
    8590                 /*Enforce D < max_damage and D > 0 */
    8591                 if(values[i]>max_damage) values[i]=max_damage;
    8592                 else if(values[i]<0) values[i]=0;
    8593         }
    8594 
    8595         /*Get all inputs and parameters*/
    8596         this->material->inputs->AddInput(new TriaInput(DamageDbarEnum,values,P1Enum));
    8597 
    8598         /*Free ressources:*/
    8599         xDelete<int>(doflist);
    8600 }
    8601 /*}}}*/
    86028057#endif
    86038058

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r16778 r16783  
    5555                /*}}}*/
    5656                /*Update virtual functions resolution: {{{*/
    57                 void  InputUpdateFromSolution(IssmDouble* solutiong);
    5857                void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
    5958                #ifdef _HAVE_DAKOTA_
    … …  
    266265                void             InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
    267266                void             InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){_error_("not implemented yet");};
    268                 void             InputUpdateFromSolutionMasstransport(IssmDouble* solution);
    269267                bool             IsInput(int name);
    270268                Gauss*         NewGauss(void);
    … …  
    295293                ElementMatrix* CreateJacobianStressbalanceSSA(void);
    296294                IssmDouble     GetYcoord(GaussTria* gauss);
    297                 void             InputUpdateFromSolutionStressbalanceHoriz( IssmDouble* solution);
    298                 void             InputUpdateFromSolutionStressbalanceFS( IssmDouble* solution);
    299                 void             InputUpdateFromSolutionStressbalanceSIA( IssmDouble* solution);
    300295                #endif
    301296
    … …  
    305300                ElementVector* CreatePVectorAdjointHoriz(void);
    306301                ElementVector* CreatePVectorAdjointBalancethickness(void);
    307                 void             InputUpdateFromSolutionAdjointHoriz( IssmDouble* solution);
    308302                #endif
    309303
    … …  
    326320                ElementVector* CreatePVectorL2ProjectionEPL(void);
    327321                void           CreateHydrologyWaterVelocityInput(void);
    328                 void             InputUpdateFromSolutionHydrology(IssmDouble* solution);
    329                 void             InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution);
    330                 void           InputUpdateFromSolutionHydrologyDC(IssmDouble* solution);
    331                 void             InputUpdateFromSolutionHydrologyDCInefficient(IssmDouble* solution);
    332                 void             InputUpdateFromSolutionHydrologyDCEfficient(IssmDouble* solution);
    333322                void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    334323                void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index);
    … …  
    348337                ElementVector* CreatePVectorDamageEvolution_CG(void);
    349338                void           DamageEvolutionF(IssmDouble* flist);
    350                 void             InputUpdateFromSolutionDamageEvolution(IssmDouble* solution);
    351339                #endif
    352340

  • issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h

    r16125 r16783  
    5454                void InputUpdateFromConstant(int constant, int name){/*Do nothing*/};
    5555                void InputUpdateFromConstant(bool constant, int name){_error_("Not implemented yet!");}
    56                 void InputUpdateFromSolution(IssmDouble* solution){_error_("Not implemented yet!");}
    5756                void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
    5857                /*}}}*/

  • issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp

    r16600 r16783  
    407407}
    408408/*}}}*/
    409 /*FUNCTION Pengrid::InputUpdateFromSolution{{{*/
    410 void  Pengrid::InputUpdateFromSolution(IssmDouble* solution){
    411         /*Nothing updated yet*/
    412 }
    413 /*}}}*/         
    414409
    415410/*Pengrid management:*/

  • issm/trunk-jpl/src/c/classes/Loads/Pengrid.h

    r16228 r16783  
    6464                void  InputUpdateFromConstant(int constant, int name);
    6565                void  InputUpdateFromConstant(bool constant, int name);
    66                 void  InputUpdateFromSolution(IssmDouble* solution);
    6766                void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
    6867                /*}}}*/

  • issm/trunk-jpl/src/c/classes/Loads/Penpair.h

    r16125 r16783  
    4444                void  InputUpdateFromConstant(int constant, int name);
    4545                void  InputUpdateFromConstant(bool constant, int name);
    46                 void  InputUpdateFromSolution(IssmDouble* solution){_error_("Not implemented yet!");}
    4746                void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
    4847                /*}}}*/

  • issm/trunk-jpl/src/c/classes/Loads/Riftfront.h

    r16237 r16783  
    6565                void    InputUpdateFromConstant(int constant, int name){_error_("Not implemented yet!");}
    6666                void    InputUpdateFromConstant(bool constant, int name);
    67                 void    InputUpdateFromSolution(IssmDouble* solution){_error_("Not implemented yet!");}
    68                 void  InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
     67                void    InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
    6968                /*}}}*/
    7069                /*Load virtual functions definitions: {{{*/

  • issm/trunk-jpl/src/c/classes/Materials/Matice.cpp

    r16638 r16783  
    854854}
    855855/*}}}*/
    856 /*FUNCTION Matice::InputUpdateFromSolution{{{*/
    857 void  Matice::InputUpdateFromSolution(IssmDouble* solution){
    858         /*Nothing updated yet*/
    859 }
    860 /*}}}*/
    861856/*FUNCTION Matice::InputUpdateFromIoModel{{{*/
    862857void Matice::InputUpdateFromIoModel(int index, IoModel* iomodel){

  • issm/trunk-jpl/src/c/classes/Materials/Matice.h

    r16638 r16783  
    4545                void  InputUpdateFromConstant(int constant, int name);
    4646                void  InputUpdateFromConstant(bool constant, int name);
    47                 void  InputUpdateFromSolution(IssmDouble* solution);
    4847                void  InputUpdateFromIoModel(int index, IoModel* iomodel);
    4948                /*}}}*/

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

    r16725 r16783  
    220220}
    221221/*}}}*/
    222 /*FUNCTION Matpar::InputUpdateFromSolution{{{*/
    223 void   Matpar::InputUpdateFromSolution(IssmDouble* solution){
    224         /*Nothing updated yet*/
    225 }
    226 /*}}}*/
    227222
    228223/*Matpar management: */

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.h

    r16725 r16783  
    7777                void   InputUpdateFromConstant(int constant, int name);
    7878                void   InputUpdateFromConstant(bool constant, int name);
    79                 void   InputUpdateFromSolution(IssmDouble* solution);
    8079                void   InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
    8180                /*}}}*/

  • issm/trunk-jpl/src/c/classes/Update.h

    r16125 r16783  
    2222                virtual void  InputUpdateFromConstant(int constant, int name)=0;
    2323                virtual void  InputUpdateFromConstant(bool constant, int name)=0;
    24                 virtual void  InputUpdateFromSolution(IssmDouble* solution)=0;
    2524                virtual void  InputUpdateFromIoModel(int index, IoModel* iomodel)=0;
    2625

  • issm/trunk-jpl/src/c/modules/InputUpdateFromSolutionx/InputUpdateFromSolutionx.cpp

    r16760 r16783  
    2828        for(int i=0;i<femmodel->elements->Size();i++){
    2929                Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
    30                 //analysis->InputUpdateFromSolution(solution,element);
    31                 element->InputUpdateFromSolution(solution);
     30                analysis->InputUpdateFromSolution(solution,element);
     31                //element->InputUpdateFromSolution(solution);
    3232        }
    3333        delete analysis;
Note: See TracChangeset for help on using the changeset viewer.