Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 482) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 483) @@ -1022,6 +1022,31 @@ int DataSet::RiftIsPresent(){ - return 0; //for now - + int i; + + vector::iterator object; + Penpair* penpair=NULL; + int found=0; + int mpi_found=0; + + /*go though loads, and figure out if one of the loads is a PenPair with numdof=2: */ + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if((*object)->Enum()==PenpairEnum()){ + + penpair=(Penpair*)(*object); + if (penpair->GetNumDofs()==2){ + found=1; + break; + } + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + return found; } @@ -1038,7 +1063,24 @@ int DataSet::MeltingIsPresent(){ - - return 0; //for now - + int found=0; + int mpi_found=0; + + vector::iterator object; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if((*object)->Enum()==PengridEnum()){ + found=1; + break; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&found,&mpi_found,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&mpi_found,1,MPI_INT,0,MPI_COMM_WORLD); + found=mpi_found; + #endif + + return found; } @@ -1047,6 +1089,34 @@ void DataSet::MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type,int sub_analysis_type){ - throw ErrorException(__FUNCT__," not implemented yet!"); - + /* generic object pointer: */ + vector::iterator object; + Pengrid* pengrid=NULL; + + int unstable; + int num_unstable_constraints=0; + int sum_num_unstable_constraints=0; + + num_unstable_constraints=0; + + /*Enforce constraints: */ + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if((*object)->Enum()==PengridEnum()){ + + pengrid=(Pengrid*)(*object); + pengrid->PenaltyConstrain(&unstable,inputs,analysis_type,sub_analysis_type); + + num_unstable_constraints+=unstable; + } + } + + #ifdef _PARALLEL_ + MPI_Reduce (&num_unstable_constraints,&sum_num_unstable_constraints,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD ); + MPI_Bcast(&sum_num_unstable_constraints,1,MPI_INT,0,MPI_COMM_WORLD); + num_unstable_constraints=sum_num_unstable_constraints; + #endif + + /*Assign output pointers: */ + *pnum_unstable_constraints=num_unstable_constraints; } Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 482) +++ /issm/trunk/src/c/Makefile.am (revision 483) @@ -188,4 +188,12 @@ ./ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp\ ./ModelProcessorx/Control/CreateParametersControl.cpp\ + ./ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp\ + ./ModelProcessorx/Thermal/CreateConstraintsThermal.cpp\ + ./ModelProcessorx/Thermal/CreateLoadsThermal.cpp\ + ./ModelProcessorx/Thermal/CreateParametersThermal.cpp\ + ./ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp\ + ./ModelProcessorx/Melting/CreateConstraintsMelting.cpp\ + ./ModelProcessorx/Melting/CreateLoadsMelting.cpp\ + ./ModelProcessorx/Melting/CreateParametersMelting.cpp\ ./Dofx/Dofx.h\ ./Dofx/Dofx.cpp\ @@ -434,4 +442,12 @@ ./ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp\ ./ModelProcessorx/Control/CreateParametersControl.cpp\ + ./ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp\ + ./ModelProcessorx/Thermal/CreateConstraintsThermal.cpp\ + ./ModelProcessorx/Thermal/CreateLoadsThermal.cpp\ + ./ModelProcessorx/Thermal/CreateParametersThermal.cpp\ + ./ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp\ + ./ModelProcessorx/Melting/CreateConstraintsMelting.cpp\ + ./ModelProcessorx/Melting/CreateLoadsMelting.cpp\ + ./ModelProcessorx/Melting/CreateParametersMelting.cpp\ ./Dofx/Dofx.h\ ./Dofx/Dofx.cpp\ Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 482) +++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 483) @@ -73,18 +73,14 @@ else if (strcmp(model->analysis_type,"thermal")==0){ - if ((strcmp(model->sub_analysis_type,"steady")==0) || (strcmp(model->sub_analysis_type,"transient")==0)){ - - //CreateElementsNodesAndMaterialsThermal(pelements,pnodes,pmaterials, model,model_handle); - //CreateConstraintsThermal(pconstraints,model,model_handle); - //CreateLoadsThermal(ploads,model,model_handle); + CreateElementsNodesAndMaterialsThermal(pelements,pnodes,pmaterials, model,model_handle); + CreateConstraintsThermal(pconstraints,model,model_handle); + CreateLoadsThermal(ploads,model,model_handle); - } - else if (strcmp(model->sub_analysis_type,"melting")==0){ + } + else if (strcmp(model->analysis_type,"melting")==0){ - //CreateElementsNodesAndMaterialsMelting(pelements,pnodes,pmaterials, model,model_handle); - //CreateConstraintsMelting(pconstraints,model,model_handle); - //CreateLoadsMelting(ploads,model,model_handle); - - } + CreateElementsNodesAndMaterialsMelting(pelements,pnodes,pmaterials, model,model_handle); + CreateConstraintsMelting(pconstraints,model,model_handle); + CreateLoadsMelting(ploads,model,model_handle); } else if (strcmp(model->analysis_type,"prognostic")==0){ Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 482) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 483) @@ -62,4 +62,5 @@ double tria_melting[3]; double tria_accumulation[3]; + double tria_geothermalflux[3]; int tria_friction_type; double tria_p; @@ -293,5 +294,5 @@ /*Create tria element using its constructor:*/ - tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot); + tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot); /*Add tria element to elements dataset: */ Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp (revision 482) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp (revision 483) @@ -43,4 +43,5 @@ /*pengrid intermediary data: */ int pengrid_id; + int pengrid_mparid; int pengrid_node_id; int pengrid_dof; @@ -141,6 +142,7 @@ pengrid_dof=1; pengrid_penalty_offset=model->penalty_offset; + pengrid_mparid=model->numberofelements+1;//refers to the corresponding parmat property card - pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate); + pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_mparid,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate); loads->AddObject(pengrid); Index: /issm/trunk/src/c/ModelProcessorx/Melting/CreateConstraintsMelting.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Melting/CreateConstraintsMelting.cpp (revision 483) +++ /issm/trunk/src/c/ModelProcessorx/Melting/CreateConstraintsMelting.cpp (revision 483) @@ -0,0 +1,33 @@ +/* + * CreateConstraintsMelting.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateConstraintsMelting" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../Model.h" + + +void CreateConstraintsMelting(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){ + + + int i; + int count; + + DataSet* constraints = NULL; + + /*Create constraints: */ + constraints = new DataSet(ConstraintsEnum()); + + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + constraints->Presort(); + + /*Assign output pointer: */ + *pconstraints=constraints; +} Index: /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp (revision 483) +++ /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp (revision 483) @@ -0,0 +1,497 @@ +/* + * CreateElementsNodesAndMaterialsMelting.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateElementsNodesAndMaterialsMelting" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../Model.h" + + +void CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){ + + + /*output: int* epart, int* my_grids, double* my_bordergrids*/ + + + int i,j,k,n; + extern int my_rank; + extern int num_procs; + + /*DataSets: */ + DataSet* elements = NULL; + DataSet* nodes = NULL; + DataSet* materials = NULL; + + /*Objects: */ + Node* node = NULL; + Penta* penta = NULL; + Matice* matice = NULL; + Matpar* matpar = NULL; + + int analysis_type; + int sub_analysis_type; + + /*output: */ + int* epart=NULL; //element partitioning. + int* npart=NULL; //node partitioning. + int* my_grids=NULL; + double* my_bordergrids=NULL; + + + /*intermediary: */ + int elements_width; //size of elements + double B_avg; + + /*matice constructor input: */ + int matice_mid; + double matice_B; + double matice_n; + + /*penta constructor input: */ + + int penta_id; + int penta_mid; + int penta_mparid; + int penta_g[6]; + double penta_h[6]; + double penta_s[6]; + double penta_b[6]; + double penta_k[6]; + int penta_friction_type; + double penta_p; + double penta_q; + int penta_shelf; + int penta_onbed; + int penta_onsurface; + double penta_meanvel;/*!scaling ratio for velocities*/ + double penta_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ + int penta_collapse; + double penta_melting[6]; + double penta_accumulation[6]; + double penta_geothermalflux[6]; + int penta_artdiff; + int penta_thermal_steadystate; + double penta_viscosity_overshoot; + double penta_stokesreconditioning; + + /*matpar constructor input: */ + int matpar_mid; + double matpar_rho_ice; + double matpar_rho_water; + double matpar_heatcapacity; + double matpar_thermalconductivity; + double matpar_latentheat; + double matpar_beta; + double matpar_meltingpoint; + double matpar_mixed_layer_capacity; + double matpar_thermal_exchange_velocity; + double matpar_g; + + /* node constructor input: */ + int node_id; + int node_partitionborder=0; + double node_x[3]; + int node_onbed; + int node_onsurface; + int node_upper_node_id; + int node_numdofs; + + + #ifdef _PARALLEL_ + /*Metis partitioning: */ + int range; + Vec gridborder=NULL; + int my_numgrids; + int* all_numgrids=NULL; + int gridcount; + int count; + #endif + int first_grid_index; + + /*Rifts:*/ + int* riftsnumpenaltypairs; + double** riftspenaltypairs; + int* riftsfill; + double* riftsfriction; + double* riftpenaltypairs=NULL; + int el1,el2; + + /*First create the elements, nodes and material properties: */ + elements = new DataSet(ElementsEnum()); + nodes = new DataSet(NodesEnum()); + materials = new DataSet(MaterialsEnum()); + + + /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ + if (!model->ismacayealpattyn)goto cleanup_and_return; + + /*Get analysis_type: */ + analysis_type=AnalysisTypeAsEnum(model->analysis_type); + sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type); + + + /*Width of elements: */ + if(strcmp(model->meshtype,"2d")==0)throw ErrorException(__FUNCT__," error message: 2d temperature computations not supported yet!"); + elements_width=6; //penta elements + + + #ifdef _PARALLEL_ + /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ + ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat"); + + + MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs); + + /*Free elements and elements2d: */ + xfree((void**)&model->elements2d); + + + /*Deal with rifts, they have to be included into one partition only, not several: */ + FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts); + + for(i=0;inumrifts;i++){ + riftpenaltypairs=model->riftspenaltypairs[i]; + for(j=0;jriftsnumpenaltypairs[i];j++){ + el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing + el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing + epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; + } + } + /*Free rifts: */ + xfree((void**)&riftsnumpenaltypairs); + for(i=0;inumrifts;i++){ + double* temp=riftspenaltypairs[i]; + xfree((void**)&temp); + } + xfree((void**)&riftspenaltypairs); + xfree((void**)&riftsfill); + xfree((void**)&riftsfriction); + + /*Used later on: */ + my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int)); + #endif + + + + /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ + + /*Fetch data needed: */ + ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); + ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat"); + ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat"); + ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat"); + ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat"); + ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat"); + ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat"); + ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat"); + ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat"); + ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat"); + ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat"); + ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat"); + ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat"); + ModelFetchData((void**)&model->accumulation,NULL,NULL,model_handle,"accumulation","Matrix","Mat"); + ModelFetchData((void**)&model->melting,NULL,NULL,model_handle,"melting","Matrix","Mat"); + + for (i=0;inumberofelements;i++){ + #ifdef _PARALLEL_ + /*We are using our element partition to decide which elements will be created on this node: */ + if(my_rank==epart[i]){ + #endif + + + /*name and id: */ + penta_id=i+1; //matlab indexing. + penta_mid=i+1; //refers to the corresponding material property card + penta_mparid=model->numberofelements+1;//refers to the corresponding parmat property card + + /*vertices,thickness,surface,bed and drag: */ + for(j=0;j<6;j++){ + penta_g[j]=(int)*(model->elements+elements_width*i+j); + penta_h[j]=*(model->thickness+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_s[j]=*(model->surface+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_b[j]=*(model->bed+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_k[j]=*(model->drag+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_melting[j]=*(model->melting+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_accumulation[j]=*(model->accumulation+ ((int)*(model->elements+elements_width*i+j)-1)); + } + + /*basal drag:*/ + penta_friction_type=(int)model->drag_type; + + penta_p=model->p[i]; + penta_q=model->q[i]; + + /*diverse: */ + penta_shelf=(int)*(model->elementoniceshelf+i); + penta_onbed=(int)*(model->elementonbed+i); + penta_onsurface=(int)*(model->elementonsurface+i); + penta_meanvel=model->meanvel; + penta_epsvel=model->epsvel; + + /*viscosity_overshoot*/ + penta_viscosity_overshoot=model->viscosity_overshoot; + penta_collapse=1; + + /*Create Penta using its constructor:*/ + penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type, + penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel, + penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff, + penta_thermal_steadystate,penta_viscosity_overshoot,penta_stokesreconditioning); + + /*Add penta element to elements dataset: */ + elements->AddObject(penta); + + + /*Deal with material:*/ + matice_mid=i+1; //same as the material id from the geom2 elements. + /*Average B over 6 element grids: */ + B_avg=0; + for(j=0;j<6;j++){ + B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1)); + } + B_avg=B_avg/6; + matice_B= B_avg; + matice_n=(double)*(model->n+i); + + /*Create matice using its constructor:*/ + matice= new Matice(matice_mid,matice_B,matice_n); + + /*Add matice element to materials dataset: */ + materials->AddObject(matice); + + #ifdef _PARALLEL_ + /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use + *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) + into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids + will hold which grids belong to this partition*/ + my_grids[(int)*(model->elements+elements_width*i+0)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+1)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+2)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+3)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+4)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+5)-1]=1; + #endif + + #ifdef _PARALLEL_ + }//if(my_rank==epart[i]) + #endif + + }//for (i=0;ielements); + xfree((void**)&model->thickness); + xfree((void**)&model->surface); + xfree((void**)&model->bed); + xfree((void**)&model->drag); + xfree((void**)&model->p); + xfree((void**)&model->q); + xfree((void**)&model->elementoniceshelf); + xfree((void**)&model->elementonbed); + xfree((void**)&model->elementonsurface); + xfree((void**)&model->elements_type); + xfree((void**)&model->n); + xfree((void**)&model->B); + + + /*Add one constant material property to materials: */ + matpar_mid=model->numberofelements+1; //put it at the end of the materials + matpar_g=model->g; + matpar_rho_ice=model->rho_ice; + matpar_rho_water=model->rho_water; + matpar_thermalconductivity=model->thermalconductivity; + matpar_heatcapacity=model->heatcapacity; + matpar_latentheat=model->latentheat; + matpar_beta=model->beta; + matpar_meltingpoint=model->meltingpoint; + matpar_mixed_layer_capacity=model->mixed_layer_capacity; + matpar_thermal_exchange_velocity=model->thermal_exchange_velocity; + + /*Create matpar object using its constructor: */ + matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity, + matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity, + matpar_thermal_exchange_velocity,matpar_g); + + /*Add to materials datset: */ + materials->AddObject(matpar); + + #ifdef _PARALLEL_ + /*From the element partitioning, we can determine which grids are on the inside of this cpu's + *element partition, and which are on its border with other nodes:*/ + gridborder=NewVec(model->numberofnodes); + + for (i=0;inumberofnodes;i++){ + if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); + } + VecAssemblyBegin(gridborder); + VecAssemblyEnd(gridborder); + + #ifdef _DEBUG_ + VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD); + #endif + + VecToMPISerial(&my_bordergrids,gridborder); + + #ifdef _DEBUG_ + if(my_rank==0){ + for (i=0;inumberofnodes;i++){ + printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]); + } + } + #endif + #endif + + /*Partition penalties in 3d: */ + if(strcmp(model->meshtype,"3d")==0){ + + /*Get penalties: */ + ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat"); + + if(model->numpenalties){ + + model->penaltypartitioning=(int*)xmalloc(model->numpenalties*sizeof(int)); + #ifdef _SERIAL_ + for(i=0;inumpenalties;i++)model->penaltypartitioning[i]=1; + #else + for(i=0;inumpenalties;i++)model->penaltypartitioning[i]=-1; + + for(i=0;inumpenalties;i++){ + first_grid_index=(int)(*(model->penalties+i*model->numlayers+0)-1); + if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition grids + /*All grids that are being penalised belong to this node's internal grid partition.:*/ + model->penaltypartitioning[i]=1; + } + if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border + model->penaltypartitioning[i]=0; + } + } + #endif + } + + /*Free penalties: */ + xfree((void**)&model->penalties); + } + + /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. + We can therefore determine which grids are internal to this node's partition + and which ones are shared with other nodes because they are on the border of this node's partition. Knowing + that, go and create the grids*/ + + /*Create nodes from x,y,z, as well as the spc values on those grids: */ + + /*First fetch data: */ + if (strcmp(model->meshtype,"3d")==0){ + ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat"); + ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat"); + } + ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat"); + ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat"); + ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat"); + ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat"); + ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat"); + + + /*Get number of dofs per node: */ + DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type); + + for (i=0;inumberofnodes;i++){ + #ifdef _PARALLEL_ + /*keep only this partition's nodes:*/ + if((my_grids[i]==1)){ + #endif + + node_id=i+1; //matlab indexing + + + + #ifdef _PARALLEL_ + if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border + node_partitionborder=1; + } + else{ + node_partitionborder=0; + } + #else + node_partitionborder=0; + #endif + + + node_x[0]=model->x[i]; + node_x[1]=model->y[i]; + node_x[2]=model->z[i]; + + + node_onbed=(int)model->gridonbed[i]; + node_onsurface=(int)model->gridonsurface[i]; + if (strcmp(model->meshtype,"3d")==0){ + if (isnan(model->uppernodes[i])){ + node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves. + } + else{ + node_upper_node_id=(int)model->uppernodes[i]; + } + } + else{ + /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/ + node_upper_node_id=node_id; + } + + /*Create node using its constructor: */ + node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id); + + /*set single point constraints.: */ + if (!model->gridonbed[i]){ + for(k=1;k<=node_numdofs;k++){ + node->FreezeDof(k); + } + } + /*Add node to nodes dataset: */ + nodes->AddObject(node); + + #ifdef _PARALLEL_ + } //if((my_grids[i]==1)) + #endif + } + + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + elements->Presort(); + nodes->Presort(); + materials->Presort(); + + /*Clean fetched data: */ + xfree((void**)&model->deadgrids); + xfree((void**)&model->x); + xfree((void**)&model->y); + xfree((void**)&model->z); + xfree((void**)&model->gridonbed); + xfree((void**)&model->gridonsurface); + xfree((void**)&model->uppernodes); + + cleanup_and_return: + + /*Assign output pointer: */ + *pelements=elements; + *pnodes=nodes; + *pmaterials=materials; + + /*Keep partitioning information into model*/ + model->epart=epart; + model->my_grids=my_grids; + model->my_bordergrids=my_bordergrids; + + /*Free ressources:*/ + #ifdef _PARALLEL_ + xfree((void**)&all_numgrids); + xfree((void**)&npart); + VecFree(&gridborder); + #endif + +} Index: /issm/trunk/src/c/ModelProcessorx/Melting/CreateLoadsMelting.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Melting/CreateLoadsMelting.cpp (revision 483) +++ /issm/trunk/src/c/ModelProcessorx/Melting/CreateLoadsMelting.cpp (revision 483) @@ -0,0 +1,87 @@ +/*! \file CreateLoadsMelting.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateLoadsMelting" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../include/macros.h" +#include "../Model.h" + + +void CreateLoadsMelting(DataSet** ploads, Model* model,ConstDataHandle model_handle){ + + int i,j,counter; + int element; + + extern int my_rank; + extern int num_procs; + + DataSet* loads = NULL; + Pengrid* pengrid = NULL; + + int segment_width; + int i1,i2,i3,i4; + int i0,range; + int grid1,grid2; + + /*pengrid intermediary data: */ + int pengrid_id; + int pengrid_node_id; + int pengrid_mparid; + int pengrid_dof; + double pengrid_penalty_offset; + int pengrid_active=0; + int pengrid_thermal_steadystate=1; + + int count; + + /*Create loads: */ + loads = new DataSet(LoadsEnum()); + count=0; + + //create penalties for grids: no grid can have a temperature over the melting point + ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat"); + + for (i=0;inumberofnodes;i++){ + #ifdef _PARALLEL_ + /*keep only this partition's nodes:*/ + if((model->my_grids[i]==1)){ + #endif + + if (model->gridonbed[i]){ + + pengrid_id=count+1; //matlab indexing + pengrid_node_id=i+1; + pengrid_dof=1; + pengrid_penalty_offset=model->penalty_offset; + pengrid_active=0; + pengrid_mparid=model->numberofelements+1;//refers to the corresponding parmat property card + + pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_mparid,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate); + + loads->AddObject(pengrid); + count++; + } + #ifdef _PARALLEL_ + } //if((model->my_grids[i]==1)) + #endif + } + xfree((void**)&model->gridonbed); + + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + loads->Presort(); + + cleanup_and_return: + + /*Assign output pointer: */ + *ploads=loads; + +} + + Index: /issm/trunk/src/c/ModelProcessorx/Melting/CreateParametersMelting.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Melting/CreateParametersMelting.cpp (revision 483) +++ /issm/trunk/src/c/ModelProcessorx/Melting/CreateParametersMelting.cpp (revision 483) @@ -0,0 +1,69 @@ +/*!\file: CreateParametersDiagnosticHoriz.cpp + * \brief driver for creating parameters dataset, for diagnostic horiz analysis. + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateParametersDiagnosticHoriz" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../Model.h" + +void CreateParametersDiagnosticHoriz(DataSet** pparameters,Model* model,ConstDataHandle model_handle){ + + Param* param = NULL; + DataSet* parameters=NULL; + int count; + int i; + + double* vx=NULL; + double* vy=NULL; + double* vz=NULL; + + /*recover parameters : */ + parameters=*pparameters; + + count=parameters->Size(); + + /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ + if (!model->ismacayealpattyn)return; + + /*Get vx and vy: */ + ModelFetchData((void**)&vx,NULL,NULL,model_handle,"vx","Matrix","Mat"); + ModelFetchData((void**)&vy,NULL,NULL,model_handle,"vy","Matrix","Mat"); + ModelFetchData((void**)&vz,NULL,NULL,model_handle,"vz","Matrix","Mat"); + + if(vx) for(i=0;inumberofnodes;i++)vx[i]=vx[i]/model->yts; + if(vy) for(i=0;inumberofnodes;i++)vy[i]=vy[i]/model->yts; + if(vz) for(i=0;inumberofnodes;i++)vz[i]=vz[i]/model->yts; + + count++; + param= new Param(count,"vx",DOUBLEVEC); + if(vx) param->SetDoubleVec(vx,model->numberofnodes); + else param->SetDoubleVec(vx,0); + parameters->AddObject(param); + + count++; + param= new Param(count,"vy",DOUBLEVEC); + if(vy) param->SetDoubleVec(vy,model->numberofnodes); + else param->SetDoubleVec(vy,0); + parameters->AddObject(param); + + count++; + param= new Param(count,"vz",DOUBLEVEC); + if(vz) param->SetDoubleVec(vz,model->numberofnodes); + else param->SetDoubleVec(vz,0); + parameters->AddObject(param); + + xfree((void**)&vx); + xfree((void**)&vy); + xfree((void**)&vz); + + /*Assign output pointer: */ + *pparameters=parameters; +} + + Index: /issm/trunk/src/c/ModelProcessorx/Model.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.h (revision 482) +++ /issm/trunk/src/c/ModelProcessorx/Model.h (revision 483) @@ -210,4 +210,16 @@ void CreateParametersControl(DataSet** pparameters,Model* model,ConstDataHandle model_handle); + /*thermal:*/ + void CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle); + void CreateConstraintsThermal(DataSet** pconstraints,Model* model,ConstDataHandle model_handle); + void CreateLoadsThermal(DataSet** ploads, Model* model, ConstDataHandle model_handle); + void CreateParametersThermal(DataSet** pparameters,Model* model,ConstDataHandle model_handle); + + /*melting:*/ + void CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle); + void CreateConstraintsMelting(DataSet** pconstraints,Model* model,ConstDataHandle model_handle); + void CreateLoadsMelting(DataSet** ploads, Model* model, ConstDataHandle model_handle); + void CreateParametersMelting(DataSet** pparameters,Model* model,ConstDataHandle model_handle); + #endif /* _MODEL_H */ Index: /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 482) +++ /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 483) @@ -58,4 +58,5 @@ double tria_melting[3]; double tria_accumulation[3]; + double tria_geothermalflux[3]; int tria_friction_type; double tria_p; @@ -226,5 +227,5 @@ /*Create tria element using its constructor:*/ - tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot); + tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot); /*Add tria element to elements dataset: */ Index: /issm/trunk/src/c/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp (revision 483) +++ /issm/trunk/src/c/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp (revision 483) @@ -0,0 +1,81 @@ +/* + * CreateConstraintsThermal.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateConstraintsThermal" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../Model.h" + + +void CreateConstraintsThermal(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){ + + + int i; + int count; + + DataSet* constraints = NULL; + Spc* spc = NULL; + + /*spc intermediary data: */ + int spc_sid; + int spc_node; + int spc_dof; + double spc_value; + + double* dirichletvalues_thermal=NULL; + double* gridondirichlet_thermal=NULL; + + /*Create constraints: */ + constraints = new DataSet(ConstraintsEnum()); + + /*Fetch data: */ + ModelFetchData((void**)&gridondirichlet_thermal,NULL,NULL,model_handle,"gridondirichlet_thermal","Matrix","Mat"); + ModelFetchData((void**)&dirichletvalues_thermal,NULL,NULL,model_handle,"dirichletvalues_thermal","Matrix","Mat"); + + count=0; + + /*Create spcs from x,y,z, as well as the spc values on those spcs: */ + for (i=0;inumberofnodes;i++){ + #ifdef _PARALLEL_ + /*keep only this partition's nodes:*/ + if((model->my_grids[i]==1)){ + #endif + + if ((int)gridondirichlet_thermal[i]){ + + /*This grid needs to be spc'd to vx_obs and vy_obs:*/ + + spc_sid=count; + spc_node=i+1; + spc_dof=1; //we enforce first translation degree of freedom, for temperature + spc_value=dirichletvalues_thermal[i]; + + spc = new Spc(spc_sid,spc_node,spc_dof,spc_value); + constraints->AddObject(spc); + count++; + } + + #ifdef _PARALLEL_ + } //if((my_grids[i]==1)) + #endif + } + + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + constraints->Presort(); + + + cleanup_and_return: + /*Free data: */ + xfree((void**)&gridondirichlet_thermal); + xfree((void**)&dirichletvalues_thermal); + + /*Assign output pointer: */ + *pconstraints=constraints; +} Index: /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp (revision 483) +++ /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp (revision 483) @@ -0,0 +1,490 @@ +/* + * CreateElementsNodesAndMaterialsThermal.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateElementsNodesAndMaterialsThermal" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../Model.h" + +#undef __FUNCT__ +#define __FUNCT__ "CreateElementsNodesAndMaterialsThermal" +void CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){ + + + /*output: int* epart, int* my_grids, double* my_bordergrids*/ + + + int i,j,k,n; + extern int my_rank; + extern int num_procs; + + /*DataSets: */ + DataSet* elements = NULL; + DataSet* nodes = NULL; + DataSet* materials = NULL; + + /*Objects: */ + Node* node = NULL; + Penta* penta = NULL; + Matice* matice = NULL; + Matpar* matpar = NULL; + + int analysis_type; + int sub_analysis_type; + + /*output: */ + int* epart=NULL; //element partitioning. + int* npart=NULL; //node partitioning. + int* my_grids=NULL; + double* my_bordergrids=NULL; + + + /*intermediary: */ + int elements_width; //size of elements + double B_avg; + + /*matice constructor input: */ + int matice_mid; + double matice_B; + double matice_n; + + /*penta constructor input: */ + + int penta_id; + int penta_mid; + int penta_mparid; + int penta_g[6]; + double penta_h[6]; + double penta_s[6]; + double penta_b[6]; + double penta_k[6]; + int penta_friction_type; + double penta_p; + double penta_q; + int penta_shelf; + int penta_onbed; + int penta_onsurface; + double penta_meanvel;/*!scaling ratio for velocities*/ + double penta_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ + int penta_collapse; + double penta_melting[6]; + double penta_accumulation[6]; + double penta_geothermalflux[6]; + int penta_artdiff; + int penta_thermal_steadystate; + double penta_viscosity_overshoot; + double penta_stokesreconditioning; + + /*matpar constructor input: */ + int matpar_mid; + double matpar_rho_ice; + double matpar_rho_water; + double matpar_heatcapacity; + double matpar_thermalconductivity; + double matpar_latentheat; + double matpar_beta; + double matpar_meltingpoint; + double matpar_mixed_layer_capacity; + double matpar_thermal_exchange_velocity; + double matpar_g; + + /* node constructor input: */ + int node_id; + int node_partitionborder=0; + double node_x[3]; + int node_onbed; + int node_onsurface; + int node_upper_node_id; + int node_numdofs; + + + #ifdef _PARALLEL_ + /*Metis partitioning: */ + int range; + Vec gridborder=NULL; + int my_numgrids; + int* all_numgrids=NULL; + int gridcount; + int count; + #endif + int first_grid_index; + + /*Rifts:*/ + int* riftsnumpenaltypairs; + double** riftspenaltypairs; + int* riftsfill; + double* riftsfriction; + double* riftpenaltypairs=NULL; + int el1,el2; + + /*First create the elements, nodes and material properties: */ + elements = new DataSet(ElementsEnum()); + nodes = new DataSet(NodesEnum()); + materials = new DataSet(MaterialsEnum()); + + + /*Get analysis_type: */ + analysis_type=AnalysisTypeAsEnum(model->analysis_type); + sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type); + + /*Width of elements: */ + if(strcmp(model->meshtype,"2d")==0){ + throw ErrorException(__FUNCT__," 2d temperature computation not supported yet!"); + } + + elements_width=6; //penta elements + + #ifdef _PARALLEL_ + /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ + ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat"); + + MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs); + + /*Free elements and elements2d: */ + xfree((void**)&model->elements2d); + + /*Deal with rifts, they have to be included into one partition only, not several: */ + FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts); + + for(i=0;inumrifts;i++){ + riftpenaltypairs=model->riftspenaltypairs[i]; + for(j=0;jriftsnumpenaltypairs[i];j++){ + el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing + el2=(int)*(riftpenaltypairs+7*j+3)-1; //matlab indexing to c indexing + epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; + } + } + /*Free rifts: */ + xfree((void**)&riftsnumpenaltypairs); + for(i=0;inumrifts;i++){ + double* temp=riftspenaltypairs[i]; + xfree((void**)&temp); + } + xfree((void**)&riftspenaltypairs); + xfree((void**)&riftsfill); + xfree((void**)&riftsfriction); + + /*Used later on: */ + my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int)); + #endif + + + /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ + + /*Fetch data needed: */ + ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); + ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat"); + ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat"); + ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat"); + ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat"); + ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat"); + ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat"); + ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat"); + ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat"); + ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat"); + ModelFetchData((void**)&model->geothermalflux,NULL,NULL,model_handle,"geothermalflux","Matrix","Mat"); + ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat"); + ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat"); + + for (i=0;inumberofelements;i++){ + #ifdef _PARALLEL_ + /*We are using our element partition to decide which elements will be created on this node: */ + if(my_rank==epart[i]){ + #endif + + + /*name and id: */ + penta_id=i+1; //matlab indexing. + penta_mid=i+1; //refers to the corresponding material property card + penta_mparid=model->numberofelements+1;//refers to the corresponding parmat property card + + /*vertices,thickness,surface,bed and drag: */ + for(j=0;j<6;j++){ + penta_g[j]=(int)*(model->elements+elements_width*i+j); + penta_h[j]=*(model->thickness+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_s[j]=*(model->surface+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_b[j]=*(model->bed+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_k[j]=*(model->drag+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_geothermalflux[j]=*(model->geothermalflux+ ((int)*(model->elements+elements_width*i+j)-1)); + } + + /*basal drag:*/ + penta_friction_type=(int)model->drag_type; + + penta_p=model->p[i]; + penta_q=model->q[i]; + + /*diverse: */ + penta_shelf=(int)*(model->elementoniceshelf+i); + penta_onbed=(int)*(model->elementonbed+i); + penta_onsurface=(int)*(model->elementonsurface+i); + penta_meanvel=model->meanvel; + penta_epsvel=model->epsvel; + penta_collapse=0; + + + /*Create Penta using its constructor:*/ + penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type, + penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel, + penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff, + penta_thermal_steadystate,penta_viscosity_overshoot,penta_stokesreconditioning); + + /*Add penta element to elements dataset: */ + elements->AddObject(penta); + + + /*Deal with material:*/ + matice_mid=i+1; //same as the material id from the geom2 elements. + /*Average B over 6 element grids: */ + B_avg=0; + for(j=0;j<6;j++){ + B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1)); + } + B_avg=B_avg/6; + matice_B= B_avg; + matice_n=(double)*(model->n+i); + + /*Create matice using its constructor:*/ + matice= new Matice(matice_mid,matice_B,matice_n); + + /*Add matice element to materials dataset: */ + materials->AddObject(matice); + + #ifdef _PARALLEL_ + /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use + *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) + into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids + will hold which grids belong to this partition*/ + my_grids[(int)*(model->elements+elements_width*i+0)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+1)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+2)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+3)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+4)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+5)-1]=1; + #endif + + #ifdef _PARALLEL_ + }//if(my_rank==epart[i]) + #endif + + }//for (i=0;ielements); + xfree((void**)&model->thickness); + xfree((void**)&model->surface); + xfree((void**)&model->bed); + xfree((void**)&model->drag); + xfree((void**)&model->p); + xfree((void**)&model->q); + xfree((void**)&model->elementoniceshelf); + xfree((void**)&model->elementonbed); + xfree((void**)&model->elementonsurface); + xfree((void**)&model->geothermalflux); + xfree((void**)&model->n); + xfree((void**)&model->B); + + /*Add one constant material property to materials: */ + matpar_mid=model->numberofelements+1; //put it at the end of the materials + matpar_g=model->g; + matpar_rho_ice=model->rho_ice; + matpar_rho_water=model->rho_water; + matpar_thermalconductivity=model->thermalconductivity; + matpar_heatcapacity=model->heatcapacity; + matpar_latentheat=model->latentheat; + matpar_beta=model->beta; + matpar_meltingpoint=model->meltingpoint; + matpar_mixed_layer_capacity=model->mixed_layer_capacity; + matpar_thermal_exchange_velocity=model->thermal_exchange_velocity; + + /*Create matpar object using its constructor: */ + matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity, + matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity, + matpar_thermal_exchange_velocity,matpar_g); + + /*Add to materials datset: */ + materials->AddObject(matpar); + + #ifdef _PARALLEL_ + /*From the element partitioning, we can determine which grids are on the inside of this cpu's + *element partition, and which are on its border with other nodes:*/ + gridborder=NewVec(model->numberofnodes); + + for (i=0;inumberofnodes;i++){ + if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); + } + VecAssemblyBegin(gridborder); + VecAssemblyEnd(gridborder); + + #ifdef _DEBUG_ + VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD); + #endif + + VecToMPISerial(&my_bordergrids,gridborder); + + #ifdef _DEBUG_ + if(my_rank==0){ + for (i=0;inumberofnodes;i++){ + printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]); + } + } + #endif + #endif + + /*Partition penalties in 3d: */ + if(strcmp(model->meshtype,"3d")==0){ + + /*Get penalties: */ + ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat"); + + if(model->numpenalties){ + + model->penaltypartitioning=(int*)xmalloc(model->numpenalties*sizeof(int)); + #ifdef _SERIAL_ + for(i=0;inumpenalties;i++)model->penaltypartitioning[i]=1; + #else + for(i=0;inumpenalties;i++)model->penaltypartitioning[i]=-1; + + for(i=0;inumpenalties;i++){ + first_grid_index=(int)(*(model->penalties+i*model->numlayers+0)-1); + if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition grids + /*All grids that are being penalised belong to this node's internal grid partition.:*/ + model->penaltypartitioning[i]=1; + } + if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border + model->penaltypartitioning[i]=0; + } + } + #endif + } + + /*Free penalties: */ + xfree((void**)&model->penalties); + } + + /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. + We can therefore determine which grids are internal to this node's partition + and which ones are shared with other nodes because they are on the border of this node's partition. Knowing + that, go and create the grids*/ + + /*Create nodes from x,y,z, as well as the spc values on those grids: */ + + /*First fetch data: */ + if (strcmp(model->meshtype,"3d")==0){ + ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat"); + ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat"); + } + ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat"); + ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat"); + ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat"); + ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat"); + ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat"); + + + /*Get number of dofs per node: */ + DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type); + + for (i=0;inumberofnodes;i++){ + #ifdef _PARALLEL_ + /*keep only this partition's nodes:*/ + if((my_grids[i]==1)){ + #endif + + node_id=i+1; //matlab indexing + + + + #ifdef _PARALLEL_ + if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border + node_partitionborder=1; + } + else{ + node_partitionborder=0; + } + #else + node_partitionborder=0; + #endif + + + node_x[0]=model->x[i]; + node_x[1]=model->y[i]; + node_x[2]=model->z[i]; + + + node_onbed=(int)model->gridonbed[i]; + node_onsurface=(int)model->gridonsurface[i]; + if (strcmp(model->meshtype,"3d")==0){ + if (isnan(model->uppernodes[i])){ + node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves. + } + else{ + node_upper_node_id=(int)model->uppernodes[i]; + } + } + else{ + /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/ + node_upper_node_id=node_id; + } + + /*Create node using its constructor: */ + node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_onbed,node_onsurface,node_upper_node_id); + + /*set single point constraints.: */ + if (strcmp(model->meshtype,"3d")==0){ + /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */ + if (model->deadgrids[i]){ + for(k=1;k<=node_numdofs;k++){ + node->FreezeDof(k); + } + } + } + /*Add node to nodes dataset: */ + nodes->AddObject(node); + + #ifdef _PARALLEL_ + } //if((my_grids[i]==1)) + #endif + } + + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + elements->Presort(); + nodes->Presort(); + materials->Presort(); + + /*Clean fetched data: */ + xfree((void**)&model->deadgrids); + xfree((void**)&model->x); + xfree((void**)&model->y); + xfree((void**)&model->z); + xfree((void**)&model->gridonbed); + xfree((void**)&model->gridonsurface); + xfree((void**)&model->uppernodes); + + cleanup_and_return: + + /*Assign output pointer: */ + *pelements=elements; + *pnodes=nodes; + *pmaterials=materials; + + /*Keep partitioning information into model*/ + model->epart=epart; + model->my_grids=my_grids; + model->my_bordergrids=my_bordergrids; + + /*Free ressources:*/ + #ifdef _PARALLEL_ + xfree((void**)&all_numgrids); + xfree((void**)&npart); + VecFree(&gridborder); + #endif + +} Index: /issm/trunk/src/c/ModelProcessorx/Thermal/CreateLoadsThermal.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Thermal/CreateLoadsThermal.cpp (revision 483) +++ /issm/trunk/src/c/ModelProcessorx/Thermal/CreateLoadsThermal.cpp (revision 483) @@ -0,0 +1,88 @@ +/*! \file CreateLoadsThermal.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateLoadsThermal" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../include/macros.h" +#include "../Model.h" + + +void CreateLoadsThermal(DataSet** ploads, Model* model,ConstDataHandle model_handle){ + + int i,j,counter; + int element; + + extern int my_rank; + extern int num_procs; + + DataSet* loads = NULL; + Pengrid* pengrid = NULL; + + int segment_width; + int i1,i2,i3,i4; + int i0,range; + int grid1,grid2; + + /*pengrid intermediary data: */ + int pengrid_id; + int pengrid_mparid; + int pengrid_node_id; + int pengrid_dof; + double pengrid_penalty_offset; + int pengrid_active=0; + int pengrid_thermal_steadystate=1; + + int numberofsegs_diag_stokes; + int count; + + /*Create loads: */ + loads = new DataSet(LoadsEnum()); + count=0; + + //create penalties for grids: no grid can have a temperature over the melting point + ModelFetchData((void**)&model->gridondirichlet_thermal,NULL,NULL,model_handle,"gridondirichlet_thermal","Matrix","Mat"); + + for (i=0;inumberofnodes;i++){ + #ifdef _PARALLEL_ + /*keep only this partition's nodes:*/ + if((model->my_grids[i]==1)){ + #endif + + if (!model->gridondirichlet_thermal[i]){ //No penalty applied on spc grids! + + pengrid_id=count+1; //matlab indexing + pengrid_node_id=i+1; + pengrid_mparid=model->numberofelements+1;//refers to the corresponding parmat property card + pengrid_dof=1; + pengrid_penalty_offset=model->penalty_offset; + pengrid_active=0; + + pengrid= new Pengrid(pengrid_id, pengrid_node_id,pengrid_mparid,pengrid_dof, pengrid_active, pengrid_penalty_offset,pengrid_thermal_steadystate); + + loads->AddObject(pengrid); + count++; + } + #ifdef _PARALLEL_ + } //if((model->my_grids[i]==1)) + #endif + } + xfree((void**)&model->gridondirichlet_thermal); + + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + loads->Presort(); + + cleanup_and_return: + + /*Assign output pointer: */ + *ploads=loads; + +} + + Index: /issm/trunk/src/c/ModelProcessorx/Thermal/CreateParametersThermal.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Thermal/CreateParametersThermal.cpp (revision 483) +++ /issm/trunk/src/c/ModelProcessorx/Thermal/CreateParametersThermal.cpp (revision 483) @@ -0,0 +1,69 @@ +/*!\file: CreateParametersDiagnosticHoriz.cpp + * \brief driver for creating parameters dataset, for diagnostic horiz analysis. + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateParametersDiagnosticHoriz" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../Model.h" + +void CreateParametersDiagnosticHoriz(DataSet** pparameters,Model* model,ConstDataHandle model_handle){ + + Param* param = NULL; + DataSet* parameters=NULL; + int count; + int i; + + double* vx=NULL; + double* vy=NULL; + double* vz=NULL; + + /*recover parameters : */ + parameters=*pparameters; + + count=parameters->Size(); + + /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ + if (!model->ismacayealpattyn)return; + + /*Get vx and vy: */ + ModelFetchData((void**)&vx,NULL,NULL,model_handle,"vx","Matrix","Mat"); + ModelFetchData((void**)&vy,NULL,NULL,model_handle,"vy","Matrix","Mat"); + ModelFetchData((void**)&vz,NULL,NULL,model_handle,"vz","Matrix","Mat"); + + if(vx) for(i=0;inumberofnodes;i++)vx[i]=vx[i]/model->yts; + if(vy) for(i=0;inumberofnodes;i++)vy[i]=vy[i]/model->yts; + if(vz) for(i=0;inumberofnodes;i++)vz[i]=vz[i]/model->yts; + + count++; + param= new Param(count,"vx",DOUBLEVEC); + if(vx) param->SetDoubleVec(vx,model->numberofnodes); + else param->SetDoubleVec(vx,0); + parameters->AddObject(param); + + count++; + param= new Param(count,"vy",DOUBLEVEC); + if(vy) param->SetDoubleVec(vy,model->numberofnodes); + else param->SetDoubleVec(vy,0); + parameters->AddObject(param); + + count++; + param= new Param(count,"vz",DOUBLEVEC); + if(vz) param->SetDoubleVec(vz,model->numberofnodes); + else param->SetDoubleVec(vz,0); + parameters->AddObject(param); + + xfree((void**)&vx); + xfree((void**)&vy); + xfree((void**)&vz); + + /*Assign output pointer: */ + *pparameters=parameters; +} + + Index: /issm/trunk/src/c/objects/Friction.cpp =================================================================== --- /issm/trunk/src/c/objects/Friction.cpp (revision 482) +++ /issm/trunk/src/c/objects/Friction.cpp (revision 483) @@ -72,4 +72,5 @@ double velocity_x[numgrids]; double velocity_y[numgrids]; + double velocity_z[numgrids]; double velocity_mag[numgrids]; @@ -90,7 +91,15 @@ //We need the velocity magnitude to evaluate the basal stress: - velocity_x[i]=*(friction->velocities+2*i+0); //velocities of size numgridsx2 - velocity_y[i]=*(friction->velocities+2*i+1); - velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2)); + if(strcmp(friction->element_type,"2d")){ + velocity_x[i]=*(friction->velocities+2*i+0); //velocities of size numgridsx2 + velocity_y[i]=*(friction->velocities+2*i+1); + velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2)); + } + else{ + velocity_x[i]=*(friction->velocities+3*i+0); //velocities of size numgridsx3 + velocity_y[i]=*(friction->velocities+3*i+1); + velocity_z[i]=*(friction->velocities+3*i+2); + velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2)+pow(velocity_z[i],2)); + } alpha2[i]=pow(friction->K[i],2)*pow(Neff[i],r)*pow(velocity_mag[i],(s-1)); Index: /issm/trunk/src/c/objects/Matpar.cpp =================================================================== --- /issm/trunk/src/c/objects/Matpar.cpp (revision 482) +++ /issm/trunk/src/c/objects/Matpar.cpp (revision 483) @@ -172,2 +172,29 @@ return new Matpar(*this); } + +double Matpar::GetMixedLayerCapacity(){ + return mixed_layer_capacity; +} + +double Matpar::GetThermalExchangeVelocity(){ + return thermal_exchange_velocity; +} + +double Matpar::GetHeatCapacity(){ + return heatcapacity; +} + +double Matpar::GetThermalConductivity(){ + return thermalconductivity; +} + +double Matpar::GetLatentHeat(){ + return latentheat; +} +double Matpar::GetBeta(){ + return beta; +} +double Matpar::GetMeltingPoint(){ + return meltingpoint; +} + Index: /issm/trunk/src/c/objects/Matpar.h =================================================================== --- /issm/trunk/src/c/objects/Matpar.h (revision 482) +++ /issm/trunk/src/c/objects/Matpar.h (revision 483) @@ -44,4 +44,11 @@ double GetRhoIce(); double GetRhoWater(); + double GetMixedLayerCapacity(); + double GetThermalExchangeVelocity(); + double GetHeatCapacity(); + double GetThermalConductivity(); + double GetLatentHeat(); + double GetBeta(); + double GetMeltingPoint(); Object* copy(); Index: /issm/trunk/src/c/objects/Pengrid.cpp =================================================================== --- /issm/trunk/src/c/objects/Pengrid.cpp (revision 482) +++ /issm/trunk/src/c/objects/Pengrid.cpp (revision 483) @@ -22,7 +22,8 @@ } -Pengrid::Pengrid(int pengrid_id, int pengrid_node_id,int pengrid_dof, int pengrid_active, double pengrid_penalty_offset,int pengrid_thermal_steadystate){ +Pengrid::Pengrid(int pengrid_id, int pengrid_mparid,int pengrid_node_id,int pengrid_dof, int pengrid_active, double pengrid_penalty_offset,int pengrid_thermal_steadystate){ id=pengrid_id; + mparid=pengrid_mparid; dof=pengrid_dof; active=pengrid_active; @@ -33,4 +34,6 @@ node_offset=UNDEF; node=NULL; + matpar=NULL; + matpar_offset=UNDEF; return; @@ -45,4 +48,5 @@ printf("Pengrid:\n"); printf(" id: %i\n",id); + printf(" mparid: %i\n",mparid); printf(" dof: %i\n",dof); printf(" active: %i\n",active); @@ -51,6 +55,8 @@ printf(" node_id: [%i]\n",node_id); printf(" node_offset: [%i]\n",node_offset); + printf(" matpar_offset=%i\n",matpar_offset); if(node)node->Echo(); + if(matpar)matpar->Echo(); return; } @@ -72,4 +78,5 @@ /*marshall Pengrid data: */ memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id); + memcpy(marshalled_dataset,&mparid,sizeof(mparid));marshalled_dataset+=sizeof(mparid); memcpy(marshalled_dataset,&dof,sizeof(dof));marshalled_dataset+=sizeof(dof); memcpy(marshalled_dataset,&active,sizeof(active));marshalled_dataset+=sizeof(active); @@ -78,4 +85,6 @@ memcpy(marshalled_dataset,&node_id,sizeof(node_id));marshalled_dataset+=sizeof(node_id); memcpy(marshalled_dataset,&node_offset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset); + memcpy(marshalled_dataset,&matpar,sizeof(matpar));marshalled_dataset+=sizeof(matpar); + memcpy(marshalled_dataset,&matpar_offset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset); *pmarshalled_dataset=marshalled_dataset; @@ -86,4 +95,5 @@ return sizeof(id)+ + sizeof(mparid)+ sizeof(dof)+ sizeof(active)+ @@ -92,4 +102,6 @@ sizeof(node_id)+ sizeof(node_offset)+ + +sizeof(matpar) + +sizeof(matpar_offset)+ sizeof(int); //sizeof(int) for enum type } @@ -117,6 +129,10 @@ memcpy(&node_id,marshalled_dataset,sizeof(node_id));marshalled_dataset+=sizeof(node_id); memcpy(&node_offset,marshalled_dataset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset); + memcpy(&matpar,marshalled_dataset,sizeof(matpar));marshalled_dataset+=sizeof(matpar); + memcpy(&matpar_offset,marshalled_dataset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset); + node=NULL; + matpar=NULL; /*return: */ @@ -143,11 +159,13 @@ DataSet* nodesin=NULL; + DataSet* materialsin=NULL; /*Recover pointers :*/ nodesin=(DataSet*)pnodesin; + materialsin=(DataSet*)pmaterialsin; /*Link this load with its nodes: */ ResolvePointers((Object**)&node,&node_id,&node_offset,1,nodesin); - + ResolvePointers((Object**)&matpar,&mparid,&matpar_offset,1,materialsin); } @@ -184,4 +202,13 @@ PenaltyCreateKMatrixDiagnosticStokes( Kgg,inputs,kmax,analysis_type,sub_analysis_type); + } + else if (analysis_type==ThermalAnalysisEnum()){ + + PenaltyCreateKMatrixThermal( Kgg,inputs,kmax,analysis_type,sub_analysis_type); + + } + else if (analysis_type==MeltingAnalysisEnum()){ + + PenaltyCreateKMatrixMelting( Kgg,inputs,kmax,analysis_type,sub_analysis_type); } @@ -202,6 +229,6 @@ int numberofdofspernode; - int dofs1=0; - int dofs2=1; + int dofs1[1]={0}; + int dofs2[1]={1}; double slope[2]; int found=0; @@ -217,7 +244,7 @@ /*recover slope: */ - found=inputs->Recover("bedslopex",&slope[0],1,&dofs1,numgrids,(void**)&node); + found=inputs->Recover("bedslopex",&slope[0],1,dofs1,numgrids,(void**)&node); if(!found)throw ErrorException(__FUNCT__," bedslopex needed in inputs!"); - found=inputs->Recover("bedslopey",&slope[1],1,&dofs2,numgrids,(void**)&node); + found=inputs->Recover("bedslopey",&slope[1],1,dofs2,numgrids,(void**)&node); if(!found)throw ErrorException(__FUNCT__," bedslopey needed in inputs!"); @@ -232,9 +259,95 @@ #undef __FUNCT__ +#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixThermal" +void Pengrid::PenaltyCreateKMatrixThermal(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){ + + int found=0; + + const int numgrids=1; + const int NDOF1=1; + const int numdof=numgrids*NDOF1; + double Ke[numdof][numdof]; + int doflist[numdof]; + int numberofdofspernode; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + + if(!active)return; + + /*Get dof list: */ + GetDofList(&doflist[0],&numberofdofspernode); + + Ke[0][0]=kmax*pow(10,penalty_offset); + + /*Add Ke to global matrix Kgg: */ + MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES); +} + +#undef __FUNCT__ +#define __FUNCT__ "Pengrid::PenaltyCreateKMatrixMelting" +void Pengrid::PenaltyCreateKMatrixMelting(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){ + + int found=0; + const int numgrids=1; + const int NDOF1=1; + const int numdof=numgrids*NDOF1; + double Ke[numdof][numdof]={0.0}; + int dofs1[1]={0}; + int doflist[numdof]; + int numberofdofspernode; + double meltingpoint; + + double pressure; + double temperature; + double beta,t_pmp; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node); + if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!"); + + found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node); + if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!"); + + /*Get dof list: */ + GetDofList(&doflist[0],&numberofdofspernode); + + //Compute pressure melting point + meltingpoint=matpar->GetMeltingPoint(); + beta=matpar->GetBeta(); + t_pmp=meltingpoint-beta*pressure; + + //Add penalty load + if (temperatureRecover("pressure",&pressure,1,dofs1,numgrids,(void**)&node); + if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!"); + + //Compute pressure melting point + meltingpoint=matpar->GetMeltingPoint(); + beta=matpar->GetBeta(); + t_pmp=meltingpoint-beta*pressure; + + //Add penalty load + P_terms[0]=kmax*pow(10,penalty_offset)*t_pmp; + + /*Add P_terms to global vector pg: */ + VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES); +} + +void Pengrid::PenaltyCreatePVectorMelting(Vec pg, void* vinputs, double kmax,int analysis_type,int sub_analysis_type){ + + const int numgrids=1; + const int NDOF1=1; + const int numdof=numgrids*NDOF1; + int doflist[numdof]; + double P_terms[numdof]={0.0}; + int numberofdofspernode; + int found=0; + int dofs1[1]={0}; + double pressure; + double temperature; + double melting_offset; + double meltingpoint; + double beta, heatcapacity; + double latentheat; + double t_pmp; + double dt; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + //First recover pressure,melting offset and temperature vectors + found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node); + if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!"); + + found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node); + if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!"); + + found=inputs->Recover("melting_offset",&melting_offset); + if(!found)throw ErrorException(__FUNCT__," could not find melting_offset in inputs!"); + + found=inputs->Recover("dt",&dt); + if((!found) && (sub_analysis_type==TransientAnalysisEnum()))throw ErrorException(__FUNCT__," could not find dt in inputs!"); + + + meltingpoint=matpar->GetMeltingPoint(); + beta=matpar->GetBeta(); + heatcapacity=matpar->GetHeatCapacity(); + latentheat=matpar->GetLatentHeat(); + + //Compute pressure melting point + t_pmp=meltingpoint-beta*pressure; + + //Add penalty load + //This time, the penalty must have the same value as the one used for the thermal computation + //so that the corresponding melting can be computed correctly + //In the thermal computation, we used kmax=melting_offset, and the same penalty_offset + if (temperatureRecover("pressure",&pressure,1,dofs1,numgrids,(void**)&node); + if(!found)throw ErrorException(__FUNCT__," could not find pressure in inputs!"); + + found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node); + if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!"); + + + //Compute pressure melting point + meltingpoint=matpar->GetMeltingPoint(); + beta=matpar->GetBeta(); + + t_pmp=meltingpoint-beta*pressure; + + //Figure out if temperature is over melting_point, in which case, this penalty needs to be activated. + + if (temperature>t_pmp){ + new_active=1; + } + else{ + new_active=0; + } + + + //Figure out stability of this penalty + if (active==new_active){ + unstable=0; + } + else{ + unstable=1; + } + + //Set penalty flag + active=new_active; + + //*Assign output pointers:*/ + *punstable=unstable; +} Index: /issm/trunk/src/c/objects/Pengrid.h =================================================================== --- /issm/trunk/src/c/objects/Pengrid.h (revision 482) +++ /issm/trunk/src/c/objects/Pengrid.h (revision 483) @@ -8,4 +8,5 @@ #include "./Node.h" #include "./Element.h" +#include "./Matpar.h" class Element; @@ -25,8 +26,12 @@ int node_offset; + int mparid; + Matpar* matpar; + int matpar_offset; + public: Pengrid(); - Pengrid(int id, int node_id,int dof, int active, double penalty_offset,int thermal_steadystate); + Pengrid(int id, int node_id,int mparid,int dof, int active, double penalty_offset,int thermal_steadystate); ~Pengrid(); @@ -49,5 +54,10 @@ void GetDofList(int* doflist,int* pnumberofdofspernode); Object* copy(); - + void PenaltyCreateKMatrixThermal(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type); + void PenaltyCreateKMatrixMelting(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type); + void MatparConfiguration(Matpar* matpar,int matpar_offset); + void PenaltyCreatePVectorThermal(Vec pg, void* inputs, double kmax,int analysis_type,int sub_analysis_type); + void PenaltyCreatePVectorMelting(Vec pg, void* inputs, double kmax,int analysis_type,int sub_analysis_type); + void PenaltyConstrain(int* punstable,void* inputs,int analysis_type,int sub_analysis_type); }; Index: /issm/trunk/src/c/objects/Penpair.cpp =================================================================== --- /issm/trunk/src/c/objects/Penpair.cpp (revision 482) +++ /issm/trunk/src/c/objects/Penpair.cpp (revision 483) @@ -260,3 +260,7 @@ return new Penpair(*this); } - + +int Penpair::GetNumDofs(){ + return numdofs; +} + Index: /issm/trunk/src/c/objects/Penpair.h =================================================================== --- /issm/trunk/src/c/objects/Penpair.h (revision 482) +++ /issm/trunk/src/c/objects/Penpair.h (revision 483) @@ -58,4 +58,5 @@ void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type,int sub_analysis_type); + int GetNumDofs(); Object* copy(); }; Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 482) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 483) @@ -313,4 +313,12 @@ CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type); } + else if (analysis_type==ThermalAnalysisEnum()){ + + CreateKMatrixThermal( Kgg,inputs,analysis_type,sub_analysis_type); + } + else if (analysis_type==MeltingAnalysisEnum()){ + + CreateKMatrixMelting( Kgg,inputs,analysis_type,sub_analysis_type); + } else{ throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet")); @@ -834,5 +842,4 @@ /*Compute thickness: */ - GetParameterValue(&thickness, &h[0],gauss_coord); /*Compute strain rate: */ @@ -1014,4 +1021,12 @@ CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type); } + else if (analysis_type==ThermalAnalysisEnum()){ + + CreatePVectorThermal( pg,inputs,analysis_type,sub_analysis_type); + } + else if (analysis_type==MeltingAnalysisEnum()){ + + CreatePVectorMelting( pg,inputs,analysis_type,sub_analysis_type); + } else{ throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet")); @@ -1208,4 +1223,5 @@ double tria_melting[3]; double tria_accumulation[3]; + double tria_geothermalflux[3]; /*configuration: */ @@ -1238,4 +1254,8 @@ tria_accumulation[2]=accumulation[g2]; + tria_geothermalflux[0]=geothermalflux[g0]; + tria_geothermalflux[1]=geothermalflux[g1]; + tria_geothermalflux[2]=geothermalflux[g2]; + tria_node_ids[0]=node_ids[g0]; tria_node_ids[1]=node_ids[g1]; @@ -1250,5 +1270,5 @@ tria_node_offsets[2]=node_offsets[g2]; - tria= new Tria(id,mid,mparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, tria_melting, tria_accumulation, friction_type,p,q,shelf,meanvel,epsvel,viscosity_overshoot); + tria= new Tria(id,mid,mparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, tria_melting, tria_accumulation, tria_geothermalflux,friction_type,p,q,shelf,meanvel,epsvel,viscosity_overshoot); tria->NodeConfiguration(tria_node_ids,tria_nodes,tria_node_offsets); @@ -3128,2 +3148,594 @@ } + +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreateKMatrixThermal" +void Penta::CreateKMatrixThermal(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ + + /* local declarations */ + int i,j; + int found=0; + + /* node data: */ + const int numgrids=6; + const int NDOF1=1; + const int numdof=NDOF1*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdof]; + int numberofdofspernode; + + /* gaussian points: */ + int num_area_gauss,igarea,igvert; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* vert_gauss_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + double gauss_weight,area_gauss_weight,vert_gauss_weight; + double gauss_coord[4]; + + int area_order=5; + int num_vert_gauss=5; + + int dofs[3]={0,1,2}; + double dt; + int dt_exists; + + double vxvyvz_list[numgrids][3]; + double vx_list[numgrids]; + int vx_list_exists; + double u; + double vy_list[numgrids]; + int vy_list_exists; + double v; + double vz_list[numgrids]; + int vz_list_exists; + double w; + + /* element parameters: */ + int onbed; + int shelf; + int steady_state; + + /*matrices: */ + double K_terms[numdof][numdof]={0.0}; + double Ke_gaussian_conduct[numdof][numdof]; + double Ke_gaussian_advec[numdof][numdof]; + double Ke_gaussian_transient[numdof][numdof]; + double B[3][numdof]; + double Bprime[3][numdof]; + double B_conduct[3][numdof]; + double B_advec[3][numdof]; + double Bprime_advec[3][numdof]; + double L[numdof]; + double D_scalar; + double D[3][3]; + double l1l2l3[3]; + double tl1l2l3D[3]; + double tBD[3][numdof]; + double tBD_conduct[3][numdof]; + double tBD_advec[3][numdof]; + double tLD[numdof]; + + double Jdet; + + /*Material properties: */ + double gravity,rho_ice,rho_water; + double heatcapacity,thermalconductivity; + double mixed_layer_capacity,thermal_exchange_velocity; + + /*Collapsed formulation: */ + Tria* tria=NULL; + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + // /*recovre material parameters: */ + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + heatcapacity=matpar->GetHeatCapacity(); + thermalconductivity=matpar->GetThermalConductivity(); + + + /*recover extra inputs from users, dt and velocity: */ + found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); + found=inputs->Recover("dt",&dt); + if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); + + for(i=0;iCreateKMatrixThermal(Kgg,inputs, analysis_type,sub_analysis_type); + delete tria; + return; + } + + /*Add Ke_gg to global matrix Kgg: */ + MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES); + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&area_gauss_weights); + xfree((void**)&vert_gauss_weights); + xfree((void**)&vert_gauss_coord); + +} + +#undef __FUNCT__ +#define __FUNCT__ "Penta::GetB_conduct" +void Penta::GetB_conduct(double* B_conduct, double* xyz_list, double* gauss_coord){ + + /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID. + * For grid i, Bi' can be expressed in the basic coordinate system + * by: + * Bi_conduct_basic=[ dh/dx ] + * [ dh/dy ] + * [ dh/dz ] + * where h is the interpolation function for grid i. + * + * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids) + */ + + int i; + int calculationdof=3; + int numgrids=6; + int DOFPERGRID=1; + + /*Same thing in the basic coordinate system: */ + double dh1dh6_basic[calculationdof][numgrids]; + + /*Get dh1dh2dh3 in basic coordinates system : */ + GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord); + + /*Build B': */ + for (i=0;iCreateKMatrixMelting(Kgg,inputs, analysis_type,sub_analysis_type); + delete tria; + return; + } +} + +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreatePVectorThermal" +void Penta::CreatePVectorThermal( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){ + + + /*indexing: */ + int i,j; + int found=0; + + const int numgrids=6; + const int NDOF1=1; + const int numdof=numgrids*NDOF1; + int doflist[numdof]; + int numberofdofspernode; + + /*Grid data: */ + double xyz_list[numgrids][3]; + + /* gaussian points: */ + int num_area_gauss,igarea,igvert; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* vert_gauss_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + double gauss_weight,area_gauss_weight,vert_gauss_weight; + double gauss_coord[4]; + int area_order=2; + int num_vert_gauss=3; + + double dt; + double vx_list[numgrids]; + double vy_list[numgrids]; + double vz_list[numgrids]; + double vxvyvz_list[numgrids][3]; + double pressure_list[3]; + double pressure; + double temperature_list[numgrids]; + double temperature; + + /*Material properties: */ + double gravity,rho_ice,rho_water; + double mixed_layer_capacity,heatcapacity; + double beta,meltingpoint,thermal_exchange_velocity; + + /* element parameters: */ + int friction_type; + + int dofs[3]={0,1,2}; + int dofs1[1]={0}; + + /*matrices: */ + double P_terms[numdof]={0.0}; + double L[numdof]; + double l1l2l3[3]; + double alpha2_list[3]; + double basalfriction_list[3]={0.0}; + double basalfriction; + double epsilon[6]; + double epsilon_sqr[3][3]; + double epsilon_matrix[3][3]; + + double Jdet; + double viscosity; + double epsilon_eff; + double phi; + double t_pmp; + double scalar; + double scalar_def; + double scalar_ocean; + double scalar_transient; + + /*Collapsed formulation: */ + Tria* tria=NULL; + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + // /*recovre material parameters: */ + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + heatcapacity=matpar->GetHeatCapacity(); + beta=matpar->GetBeta(); + meltingpoint=matpar->GetMeltingPoint(); + + /*recover extra inputs from users, dt and velocity: */ + found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); + 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!"); + + for(i=0;iGetViscosity3d(&viscosity,&epsilon[0]); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); + + /* Get nodal functions */ + GetNodalFunctions(&L[0], gauss_coord); + + /*Build deformational heating: */ + GetPhi(&phi, &epsilon[0], viscosity); + + /*Build pe_gaussian */ + if(sub_analysis_type==SteadyAnalysisEnum()){ + scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss_weight; + } + else{ + scalar_def=dt*phi/(rho_ice*heatcapacity)*Jdet*gauss_weight; + } + + 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); +} + +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreatePVectorMelting" +void Penta::CreatePVectorMelting( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){ + return; +} + +#undef __FUNCT__ +#define __FUNCT__ "Penta::GetPhi" +void Penta::GetPhi(double* phi, double* epsilon, double viscosity){ + /*Compute deformational heating from epsilon and viscosity */ + + double epsilon_matrix[3][3]; + double epsilon_eff; + double epsilon_sqr[3][3]; + + /* Build epsilon matrix */ + epsilon_matrix[0][0]=*(epsilon+0); + epsilon_matrix[1][0]=*(epsilon+3); + epsilon_matrix[2][0]=*(epsilon+4); + epsilon_matrix[0][1]=*(epsilon+3); + epsilon_matrix[1][1]=*(epsilon+1); + epsilon_matrix[2][1]=*(epsilon+5); + epsilon_matrix[0][2]=*(epsilon+4); + epsilon_matrix[1][2]=*(epsilon+5); + epsilon_matrix[2][2]=*(epsilon+2); + + /* Effective value of epsilon_matrix */ + epsilon_sqr[0][0]=pow(epsilon_matrix[0][0],2); + epsilon_sqr[1][0]=pow(epsilon_matrix[1][0],2); + epsilon_sqr[2][0]=pow(epsilon_matrix[2][0],2); + epsilon_sqr[0][1]=pow(epsilon_matrix[0][1],2); + epsilon_sqr[1][1]=pow(epsilon_matrix[1][1],2); + epsilon_sqr[2][1]=pow(epsilon_matrix[2][1],2); + epsilon_sqr[0][2]=pow(epsilon_matrix[0][2],2); + epsilon_sqr[1][2]=pow(epsilon_matrix[1][2],2); + epsilon_sqr[2][2]=pow(epsilon_matrix[2][2],2); + + epsilon_eff=1/pow(2,1.0/2.0)*pow((epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]),1.0/2.0); + *phi=2*pow(epsilon_eff,2.0)*viscosity; +} Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 482) +++ /issm/trunk/src/c/objects/Penta.h (revision 483) @@ -130,4 +130,12 @@ void ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp); void GetNodalFunctionsStokes(double* l1l7, double* gauss_coord); + void CreateKMatrixThermal(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); + void GetB_conduct(double* B_conduct, double* xyz_list, double* gauss_coord); + void GetB_advec(double* B_advec, double* xyz_list, double* gauss_coord); + void GetBprime_advec(double* Bprime_advec, double* xyz_list, double* gauss_coord); + void CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); + void CreatePVectorThermal( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type); + void CreatePVectorMelting( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type); + void GetPhi(double* phi, double* epsilon, double viscosity); Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 482) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 483) @@ -33,5 +33,5 @@ Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_node_ids[3],double tria_h[3],double tria_s[3],double tria_b[3],double tria_k[3],double tria_melting[3], - double tria_accumulation[3],int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel, + double tria_accumulation[3],double tria_geothermalflux[3],int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel, double tria_viscosity_overshoot){ @@ -51,5 +51,5 @@ melting[i]=tria_melting[i]; accumulation[i]=tria_accumulation[i]; - } + geothermalflux[i]=tria_geothermalflux[i]; } matice=NULL; matice_offset=UNDEF; @@ -87,4 +87,5 @@ printf(" melting=[%g,%g,%g]\n",melting[0],melting[1],melting[2]); printf(" accumulation=[%g,%g,%g]\n",accumulation[0],accumulation[1],accumulation[2]); + printf(" geothermalflux=[%g,%g,%g]\n",geothermalflux[0],geothermalflux[1],geothermalflux[2]); printf(" friction_type: %i\n",friction_type); printf(" p: %g\n",p); @@ -136,4 +137,5 @@ memcpy(marshalled_dataset,&melting,sizeof(melting));marshalled_dataset+=sizeof(melting); memcpy(marshalled_dataset,&accumulation,sizeof(accumulation));marshalled_dataset+=sizeof(accumulation); + memcpy(marshalled_dataset,&geothermalflux,sizeof(geothermalflux));marshalled_dataset+=sizeof(geothermalflux); memcpy(marshalled_dataset,&friction_type,sizeof(friction_type));marshalled_dataset+=sizeof(friction_type); memcpy(marshalled_dataset,&onbed,sizeof(onbed));marshalled_dataset+=sizeof(onbed); @@ -166,4 +168,5 @@ +sizeof(melting) +sizeof(accumulation) + +sizeof(geothermalflux) +sizeof(friction_type) +sizeof(onbed) @@ -208,4 +211,5 @@ memcpy(&melting,marshalled_dataset,sizeof(melting));marshalled_dataset+=sizeof(melting); memcpy(&accumulation,marshalled_dataset,sizeof(accumulation));marshalled_dataset+=sizeof(accumulation); + memcpy(&geothermalflux,marshalled_dataset,sizeof(geothermalflux));marshalled_dataset+=sizeof(geothermalflux); memcpy(&friction_type,marshalled_dataset,sizeof(friction_type));marshalled_dataset+=sizeof(friction_type); memcpy(&onbed,marshalled_dataset,sizeof(onbed));marshalled_dataset+=sizeof(onbed); @@ -264,5 +268,5 @@ #undef __FUNCT__ -#define __FUNCT__ "Tria::CreateKMatrix" +#define __FUNCT__ "Tria::reateKMatrix" void Tria::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ @@ -305,7 +309,7 @@ /* node data: */ const int numgrids=3; - const int numdofs=2*numgrids; + const int numdof=2*numgrids; double xyz_list[numgrids][3]; - int doflist[numdofs]; + int doflist[numdof]; int numberofdofspernode; @@ -329,12 +333,12 @@ /* matrices: */ - double B[3][numdofs]; - double Bprime[3][numdofs]; + double B[3][numdof]; + double Bprime[3][numdof]; double D[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}}; // material matrix, simple scalar matrix. double D_scalar; /* local element matrices: */ - double Ke_gg[numdofs][numdofs]; //local element stiffness matrix - double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix evaluated at the gaussian point. + double Ke_gg[numdof][numdof]; //local element stiffness matrix + double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. double Jdet; @@ -360,5 +364,5 @@ /* Set Ke_gg to 0: */ - for(i=0;iRecover("drag",&k[0],1,dofs,3,(void**)nodes); inputs->Recover("melting",&melting[0],1,dofs,3,(void**)nodes); - inputs->Recover("accumulation_param",&accumulation[0],1,dofs,3,(void**)nodes); + inputs->Recover("accumulation",&accumulation[0],1,dofs,3,(void**)nodes); + inputs->Recover("geothermalflux",&geothermalflux[0],1,dofs,3,(void**)nodes); //Update material if necessary @@ -1556,8 +1560,8 @@ /* node data: */ const int numgrids=3; - const int numdofs=2*numgrids; + const int numdof=2*numgrids; const int NDOF2=2; double xyz_list[numgrids][3]; - int doflist[numdofs]; + int doflist[numdof]; int numberofdofspernode; @@ -1588,6 +1592,6 @@ /*element vector : */ - double due_g[numdofs]; - double due_g_gaussian[numdofs]; + double due_g[numdof]; + double due_g_gaussian[numdof]; /* Jacobian: */ @@ -1613,5 +1617,5 @@ /* Set due_g to 0: */ - for(i=0;iRecover("dt",&dt); + if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!"); + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + //recover material parameters + mixed_layer_capacity=matpar->GetMixedLayerCapacity(); + thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + + + GaussTria (&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gauss (nodes on the bedrock) */ + for (ig=0; igGetLatentHeat(); + heatcapacity=matpar->GetHeatCapacity(); + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Get gaussian points and weights: */ + GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gauss (nodes on the bedrock) */ + for (ig=0; igGetMixedLayerCapacity(); + thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + beta=matpar->GetBeta(); + meltingpoint=matpar->GetMeltingPoint(); + + + /*recover extra inputs from users, dt and velocity: */ + 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!"); + + /* Ice/ocean heat exchange flux on ice shelf base */ + + GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gauss 2d (nodes on the bedrock) */ + for (ig=0; igGetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + + + /*recover extra inputs from users, dt and velocity: */ + found=inputs->Recover("dt",&dt); + if((!found) && (sub_analysis_type==TransientAnalysisEnum()))throw ErrorException(__FUNCT__," could not find dt in inputs!"); + + found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find velocity in inputs!"); + + for(i=0;ielement_type=(char*)xmalloc((strlen("3d")+1)*sizeof(char)); + strcpy(friction->element_type,"3d"); + + friction->gravity=matpar->GetG(); + friction->rho_ice=matpar->GetRhoIce(); + friction->rho_water=matpar->GetRhoWater(); + friction->K=&k[0]; + friction->bed=&b[0]; + friction->thickness=&h[0]; + friction->velocities=&vxvyvz_list[0][0]; + friction->p=p; + friction->q=q; + + /*Compute alpha2_list: */ + FrictionGetAlpha2(&alpha2_list[0],friction); + + /*Erase friction object: */ + DeleteFriction(&friction); + + /* Compute basal friction */ + for(i=0;i