Index: /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp =================================================================== --- /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp (revision 502) +++ /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp (revision 503) @@ -22,4 +22,5 @@ /*intermediary: */ + elements->Configure(elements,loads,nodes,materials); loads->Configure(elements,loads,nodes,materials); Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 502) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 503) @@ -1105,4 +1105,5 @@ pengrid=(Pengrid*)(*object); + pengrid->PenaltyConstrain(&unstable,inputs,analysis_type,sub_analysis_type); Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 502) +++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 503) @@ -76,4 +76,5 @@ CreateConstraintsThermal(pconstraints,model,model_handle); CreateLoadsThermal(ploads,model,model_handle); + CreateParametersThermal(pparameters,model,model_handle); } @@ -83,4 +84,5 @@ CreateConstraintsMelting(pconstraints,model,model_handle); CreateLoadsMelting(ploads,model,model_handle); + CreateParametersMelting(pparameters,model,model_handle); } else if (strcmp(model->analysis_type,"prognostic")==0){ Index: /issm/trunk/src/c/ModelProcessorx/Melting/CreateParametersMelting.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Melting/CreateParametersMelting.cpp (revision 502) +++ /issm/trunk/src/c/ModelProcessorx/Melting/CreateParametersMelting.cpp (revision 503) @@ -1,8 +1,8 @@ -/*!\file: CreateParametersDiagnosticHoriz.cpp +/*!\file: CreateParametersMelting.cpp * \brief driver for creating parameters dataset, for diagnostic horiz analysis. */ #undef __FUNCT__ -#define __FUNCT__ "CreateParametersDiagnosticHoriz" +#define __FUNCT__ "CreateParametersMelting" #include "../../DataSet/DataSet.h" @@ -13,5 +13,5 @@ #include "../Model.h" -void CreateParametersDiagnosticHoriz(DataSet** pparameters,Model* model,ConstDataHandle model_handle){ +void CreateParametersMelting(DataSet** pparameters,Model* model,ConstDataHandle model_handle){ Param* param = NULL; Index: /issm/trunk/src/c/ModelProcessorx/Thermal/CreateParametersThermal.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Thermal/CreateParametersThermal.cpp (revision 502) +++ /issm/trunk/src/c/ModelProcessorx/Thermal/CreateParametersThermal.cpp (revision 503) @@ -1,8 +1,8 @@ -/*!\file: CreateParametersDiagnosticHoriz.cpp +/*!\file: CreateParametersThermal.cpp * \brief driver for creating parameters dataset, for diagnostic horiz analysis. */ #undef __FUNCT__ -#define __FUNCT__ "CreateParametersDiagnosticHoriz" +#define __FUNCT__ "CreateParametersThermal" #include "../../DataSet/DataSet.h" @@ -13,5 +13,5 @@ #include "../Model.h" -void CreateParametersDiagnosticHoriz(DataSet** pparameters,Model* model,ConstDataHandle model_handle){ +void CreateParametersThermal(DataSet** pparameters,Model* model,ConstDataHandle model_handle){ Param* param = NULL; @@ -23,4 +23,5 @@ double* vy=NULL; double* vz=NULL; + double* pressure=NULL; /*recover parameters : */ @@ -63,4 +64,17 @@ xfree((void**)&vz); + /*Get pressure: */ + ModelFetchData((void**)&pressure,NULL,NULL,model_handle,"pressure","Matrix","Mat"); + + count++; + param= new Param(count,"pressure",DOUBLEVEC); + if(pressure) param->SetDoubleVec(pressure,model->numberofnodes); + else param->SetDoubleVec(pressure,0); + parameters->AddObject(param); + + /*Free pressure: */ + xfree((void**)&pressure); + + /*Assign output pointer: */ *pparameters=parameters; Index: /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp =================================================================== --- /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp (revision 502) +++ /issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp (revision 503) @@ -24,4 +24,5 @@ double* vy=NULL; double* vz=NULL; + double* pressure=NULL; double* u_g=NULL; @@ -54,5 +55,5 @@ - if ( (analysis_type==ControlAnalysisEnum()) || (analysis_type==DiagnosticAnalysisEnum())){ + if ( (analysis_type==ControlAnalysisEnum()) || (analysis_type==DiagnosticAnalysisEnum()) || (analysis_type==ThermalAnalysisEnum())){ parameters->FindParam((void*)&vx,"vx"); @@ -141,4 +142,26 @@ } + if(analysis_type==ThermalAnalysisEnum()){ + + parameters->FindParam((void*)&pressure,"pressure"); + + /*Now, from pressure, build p_g, correctly partitioned: */ + p_g=(double*)xcalloc(numberofnodes,sizeof(double)); + + for(i=0;iSetDoubleVec(p_g,numberofnodes); + parameters->AddObject(param); + + /*erase old parameter: */ + param=(Param*)parameters->FindParamObject("pressure"); + parameters->DeleteObject((Object*)param); + } + xfree((void**)&partition); @@ -150,4 +173,5 @@ xfree((void**)&vx_obs); xfree((void**)&vy_obs); + xfree((void**)&pressure); xfree((void**)&control_parameter); xfree((void**)&u_g_obs); Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 502) +++ /issm/trunk/src/c/objects/Node.cpp (revision 503) @@ -444,5 +444,5 @@ ParameterInputs* inputs=NULL; Node* node=this; - int* dof={0}; + int dof[1]={0}; /*Recover parameter inputs: */ Index: /issm/trunk/src/c/objects/Pengrid.cpp =================================================================== --- /issm/trunk/src/c/objects/Pengrid.cpp (revision 502) +++ /issm/trunk/src/c/objects/Pengrid.cpp (revision 503) @@ -529,5 +529,5 @@ double meltingpoint; int new_active; - int* dofs1={0}; + int dofs1[1]={0}; int unstable=0; @@ -569,5 +569,5 @@ unstable=1; } - + //Set penalty flag active=new_active; Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 502) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 503) @@ -3539,6 +3539,4 @@ double vz_list[numgrids]; double vxvyvz_list[numgrids][3]; - double pressure_list[3]; - double pressure; double temperature_list[numgrids]; double temperature; @@ -3599,13 +3597,12 @@ if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); - found=inputs->Recover("dt",&dt); - if((!found) && (sub_analysis_type==TransientAnalysisEnum()))throw ErrorException(__FUNCT__," could not find dt in inputs!"); - - found=inputs->Recover("pressure",&pressure_list[0],1,dofs1,numgrids,(void**)nodes); - if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!"); - - found=inputs->Recover("temperature",&temperature_list[0],1,dofs1,numgrids,(void**)nodes); - if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!"); - + if(sub_analysis_type==TransientAnalysisEnum()){ + found=inputs->Recover("dt",&dt); + if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); + + found=inputs->Recover("temperature",&temperature_list[0],1,dofs1,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!"); + } + for(i=0;iCreatePVectorThermalShelf(pg,inputs, analysis_type,sub_analysis_type); delete tria; - return; } /* Geothermal flux on ice sheet base and basal friction */ if(onbed && !shelf){ - + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. tria->CreatePVectorThermalSheet(pg,inputs, analysis_type,sub_analysis_type); delete tria; - return; } cleanup_and_return: - xfree((void**)first_gauss_area_coord); - xfree((void**)second_gauss_area_coord); - xfree((void**)third_gauss_area_coord); - xfree((void**)vert_gauss_coord); - xfree((void**)area_gauss_weights); - xfree((void**)vert_gauss_weights); + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&vert_gauss_coord); + xfree((void**)&area_gauss_weights); + xfree((void**)&vert_gauss_weights); + } Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 502) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 503) @@ -2598,5 +2598,5 @@ /*recover extra inputs from users, dt: */ - inputs->Recover("dt",&dt); + found=inputs->Recover("dt",&dt); if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); Index: /issm/trunk/src/m/solutions/cielo/thermal.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/thermal.m (revision 502) +++ /issm/trunk/src/m/solutions/cielo/thermal.m (revision 503) @@ -19,5 +19,5 @@ %initialize inputs - displaystring(debug,'\n%s',['setup inputs...']); + displaystring(md.debug,'\n%s',['setup inputs...']); inputs=inputlist; inputs=add(inputs,'velocity',m_t.parameters.u_g,'doublevec',3,m_t.parameters.numberofnodes); @@ -27,13 +27,13 @@ if strcmpi(md.sub_analysis_type,'steady'), - displaystring(debug,'\n%s',['computing temperatures...']); + displaystring(md.debug,'\n%s',['computing temperatures...']); [t_g m_t.loads melting_offset]=thermal_core(m_t,inputs,'thermal','steady'); - displaystring(debug,'\n%s',['computing melting...']); + displaystring(md.debug,'\n%s',['computing melting...']); inputs=add(inputs,'melting_offset',melting_offset,'double'); inputs=add(inputs,'temperature',t_g,'doublevec',1,m_t.parameters.numberofnodes); melting_g=diagnostic_core_linear(m_m,inputs,'melting','steady'); - displaystring(debug,'\n%s',['load results...']); + displaystring(md.debug,'\n%s',['load results...']); md.temperature=t_g; md.melting=melting_g*md.yts; %from m/s to m/a @@ -53,7 +53,7 @@ for n=1:nsteps, - displaystring(debug,'\n%s%i/%i\n','time step: ',n,nsteps); + displaystring(md.debug,'\n%s%i/%i\n','time step: ',n,nsteps); - displaystring(debug,'\n%s',[' computing temperatures...']); + displaystring(md.debug,'\n%s',[' computing temperatures...']); inputs=add(inputs,'temperature',soln(n).t_g,'doublevec',1,m_t.parameters.numberofnodes); [soln(n+1).t_g m_t.loads melting_offset]=thermal_core(m_t,inputs); Index: /issm/trunk/src/m/solutions/cielo/thermal_core.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/thermal_core.m (revision 502) +++ /issm/trunk/src/m/solutions/cielo/thermal_core.m (revision 503) @@ -26,6 +26,8 @@ %system matrices - [K_gg_ , p_g_]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs); - [K_gg_ , p_g_]=SystemMatrices(K_gg,p_g,m.elements,m.nodes,loads,m.materials,m.parameters,inputs); + displaystring(m.parameters.debug,'%s',[' system matrices']); + [K_gg , p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + displaystring(m.parameters.debug,'%s',[' penalty system matrices']); + [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); %Reduce tangent matrix from g size to f size @@ -36,21 +38,25 @@ %Solve - displaystring(m.parameters.debug,'\n%s',[' condition number of stiffness matrix: ',condest(K_ff)]); + displaystring(m.parameters.debug,'%s',[' condition number of stiffness matrix: ',condest(K_ff)]); [t_f]=Solver(K_ff,p_f,[],m.parameters); %Merge back to g set + displaystring(m.parameters.debug,'%s',[' merging solution back to g set']); [t_g]= Mergesolutionfromftog( t_f, m.Gmn, m.ys, m.nodesets ); %Update inputs in datasets inputs=add(inputs,'temperature',t_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes); + displaystring(m.parameters.debug,'%s',[' update inputs']); [m.elements,m.nodes, loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,inputs); %penalty constraints - [m.loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, m.loads, m.materials,m.parameters,inputs); + displaystring(m.parameters.debug,'%s',[' penalty constraints']); + [m.loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, m.loads, m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); if ~converged, displaystring(m.parameters.debug,'\n%s%i','#unstable constraints ',num_unstable_constraints); + m.parameters.min_thermal_constraints - if num_unstable_constraints