Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 357) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 358) @@ -1153,4 +1153,20 @@ } +void DataSet::SlopeExtrude(Vec sg,double* sg_serial){ + + vector::iterator object; + Penta* penta=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if((*object)->Enum()==PentaEnum()){ + + penta=(Penta*)(*object); + penta->SlopeExtrude(sg,sg_serial); + + } + } + +} void DataSet::ComputePressure(Vec p_g){ Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 357) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 358) @@ -75,4 +75,5 @@ void Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type); void VelocityExtrude(Vec ug,double* ug_serial); + void SlopeExtrude(Vec sg,double* sg_serial); int DeleteObject(Object* object); void ComputePressure(Vec p_g); Index: /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp (revision 357) +++ /issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp (revision 358) @@ -42,9 +42,18 @@ return DiagnosticHutterAnalysisEnum(); } - else if (strcmp(analysis_type,"surface_slope_compute")==0){ - return SurfaceSlopeComputeAnalysisEnum(); + else if (strcmp(analysis_type,"slope_compute")==0){ + return SlopeComputeAnalysisEnum(); } - else if (strcmp(analysis_type,"bed_slope_compute")==0){ - return BedSlopeComputeAnalysisEnum(); + else if (strcmp(analysis_type,"surface_slope_compute_x")==0){ + return SurfaceSlopeComputeXAnalysisEnum(); + } + else if (strcmp(analysis_type,"surface_slope_compute_y")==0){ + return SurfaceSlopeComputeYAnalysisEnum(); + } + else if (strcmp(analysis_type,"bed_slope_compute_x")==0){ + return BedSlopeComputeXAnalysisEnum(); + } + else if (strcmp(analysis_type,"bed_slope_compute_y")==0){ + return BedSlopeComputeYAnalysisEnum(); } else throw ErrorException(__FUNCT__," could not recognized analysis type!"); Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 357) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 358) @@ -26,6 +26,9 @@ int DiagnosticStokesAnalysisEnum(void){ return 30; } int DiagnosticHutterAnalysisEnum(void){ return 31; } -int SurfaceSlopeComputeAnalysisEnum(void){ return 32; } -int BedSlopeComputeAnalysisEnum(void){ return 33; } +int SlopeComputeAnalysisEnum(void){ return 32; } +int SurfaceSlopeComputeXAnalysisEnum(void){ return 33; } +int SurfaceSlopeComputeYAnalysisEnum(void){ return 34; } +int BedSlopeComputeXAnalysisEnum(void){ return 35; } +int BedSlopeComputeYAnalysisEnum(void){ return 36; } /*datasets: */ Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 357) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 358) @@ -42,6 +42,9 @@ int DiagnosticStokesAnalysisEnum(void); int DiagnosticHutterAnalysisEnum(void); -int SurfaceSlopeComputeAnalysisEnum(void); -int BedSlopeComputeAnalysisEnum(void); +int SlopeComputeAnalysisEnum(void); +int SurfaceSlopeComputeXAnalysisEnum(void); +int SurfaceSlopeComputeYAnalysisEnum(void); +int BedSlopeComputeXAnalysisEnum(void); +int BedSlopeComputeYAnalysisEnum(void); Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 357) +++ /issm/trunk/src/c/Makefile.am (revision 358) @@ -179,10 +179,7 @@ ./ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp \ ./ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp\ - ./ModelProcessorx/SurfaceSlopeCompute/CreateElementsNodesAndMaterialsSurfaceSlopeCompute.cpp\ - ./ModelProcessorx/SurfaceSlopeCompute/CreateConstraintsSurfaceSlopeCompute.cpp \ - ./ModelProcessorx/SurfaceSlopeCompute/CreateLoadsSurfaceSlopeCompute.cpp\ - ./ModelProcessorx/BedSlopeCompute/CreateElementsNodesAndMaterialsBedSlopeCompute.cpp\ - ./ModelProcessorx/BedSlopeCompute/CreateConstraintsBedSlopeCompute.cpp \ - ./ModelProcessorx/BedSlopeCompute/CreateLoadsBedSlopeCompute.cpp\ + ./ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp\ + ./ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp \ + ./ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp\ ./ModelProcessorx/Control/CreateParametersControl.cpp\ ./Dofx/Dofx.h\ @@ -243,5 +240,8 @@ ./ProcessParamsx/ProcessParamsx.h\ ./VelocityExtrudex/VelocityExtrudex.cpp\ - ./VelocityExtrudex/VelocityExtrudex.h + ./VelocityExtrudex/VelocityExtrudex.h\ + ./SlopeExtrudex/SlopeExtrudex.cpp\ + ./SlopeExtrudex/SlopeExtrudex.h + @@ -420,10 +420,7 @@ ./ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp \ ./ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp\ - ./ModelProcessorx/SurfaceSlopeCompute/CreateElementsNodesAndMaterialsSurfaceSlopeCompute.cpp\ - ./ModelProcessorx/SurfaceSlopeCompute/CreateConstraintsSurfaceSlopeCompute.cpp \ - ./ModelProcessorx/SurfaceSlopeCompute/CreateLoadsSurfaceSlopeCompute.cpp\ - ./ModelProcessorx/BedSlopeCompute/CreateElementsNodesAndMaterialsBedSlopeCompute.cpp\ - ./ModelProcessorx/BedSlopeCompute/CreateConstraintsBedSlopeCompute.cpp \ - ./ModelProcessorx/BedSlopeCompute/CreateLoadsBedSlopeCompute.cpp\ + ./ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp\ + ./ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp \ + ./ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp\ ./ModelProcessorx/Control/CreateParametersControl.cpp\ ./Dofx/Dofx.h\ @@ -485,4 +482,6 @@ ./VelocityExtrudex/VelocityExtrudex.cpp\ ./VelocityExtrudex/VelocityExtrudex.h\ + ./SlopeExtrudex/SlopeExtrudex.cpp\ + ./SlopeExtrudex/SlopeExtrudex.h\ ./parallel/diagnostic_core_nonlinear.cpp\ ./parallel/CreateFemModel.cpp\ Index: /issm/trunk/src/c/ModelProcessorx/CreateConstraints.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateConstraints.cpp (revision 357) +++ /issm/trunk/src/c/ModelProcessorx/CreateConstraints.cpp (revision 358) @@ -30,9 +30,6 @@ CreateConstraintsDiagnosticHutter(pconstraints,model,model_handle); } - else if (strcmp(model->analysis_type,"surface_slope_compute")==0){ - CreateConstraintsSurfaceSlopeCompute(pconstraints,model,model_handle); - } - else if (strcmp(model->analysis_type,"bed_slope_compute")==0){ - CreateConstraintsBedSlopeCompute(pconstraints,model,model_handle); + else if (strcmp(model->analysis_type,"slope_compute")==0){ + CreateConstraintsSlopeCompute(pconstraints,model,model_handle); } /* Index: /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp (revision 357) +++ /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp (revision 358) @@ -39,12 +39,7 @@ } - else if ((strcmp(model->analysis_type,"surface_slope_compute")==0)){ + else if ((strcmp(model->analysis_type,"slope_compute")==0)){ - CreateElementsNodesAndMaterialsSurfaceSlopeCompute(pelements,pnodes,pmaterials, model,model_handle); - - } - else if ((strcmp(model->analysis_type,"bed_slope_compute")==0)){ - - CreateElementsNodesAndMaterialsBedSlopeCompute(pelements,pnodes,pmaterials, model,model_handle); + CreateElementsNodesAndMaterialsSlopeCompute(pelements,pnodes,pmaterials, model,model_handle); } Index: /issm/trunk/src/c/ModelProcessorx/CreateLoads.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateLoads.cpp (revision 357) +++ /issm/trunk/src/c/ModelProcessorx/CreateLoads.cpp (revision 358) @@ -34,11 +34,7 @@ CreateLoadsDiagnosticHutter(ploads,model,model_handle); } - else if (strcmp(model->analysis_type,"surface_slope_compute")==0){ + else if (strcmp(model->analysis_type,"slope_compute")==0){ - CreateLoadsSurfaceSlopeCompute(ploads,model,model_handle); - } - else if (strcmp(model->analysis_type,"bed_slope_compute")==0){ - - CreateLoadsBedSlopeCompute(ploads,model,model_handle); + CreateLoadsSlopeCompute(ploads,model,model_handle); } /* Index: /issm/trunk/src/c/ModelProcessorx/Model.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 357) +++ /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 358) @@ -293,5 +293,4 @@ ModelFetchData((void**)&model->ishutter,NULL,NULL,model_handle,"ishutter","Integer",NULL); ModelFetchData((void**)&model->ismacayealpattyn,NULL,NULL,model_handle,"ismacayealpattyn","Integer",NULL); - ModelFetchData((void**)&model->isstokes,NULL,NULL,model_handle,"isstokes","Integer",NULL); /*!Get drag_type, drag and p,q: */ @@ -340,5 +339,5 @@ ModelFetchData((void**)&model->mixed_layer_capacity,NULL,NULL,model_handle,"mixed_layer_capacity","Scalar",NULL); ModelFetchData((void**)&model->thermal_exchange_velocity,NULL,NULL,model_handle,"thermal_exchange_velocity","Scalar",NULL); - + /*rifts: */ ModelFetchData((void**)&model->numrifts,NULL,NULL,model_handle,"numrifts","Integer",NULL); Index: /issm/trunk/src/c/ModelProcessorx/Model.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.h (revision 357) +++ /issm/trunk/src/c/ModelProcessorx/Model.h (revision 358) @@ -200,15 +200,9 @@ void CreateParametersDiagnosticStokes(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount); - /*surface slope compute*/ - void CreateElementsNodesAndMaterialsSurfaceSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle); - void CreateConstraintsSurfaceSlopeCompute(DataSet** pconstraints,Model* model,ConstDataHandle model_handle); - void CreateLoadsSurfaceSlopeCompute(DataSet** ploads, Model* model, ConstDataHandle model_handle); - void CreateParametersSurfaceSlopeCompute(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount); - - /*bed slope compute*/ - void CreateElementsNodesAndMaterialsBedSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle); - void CreateConstraintsBedSlopeCompute(DataSet** pconstraints,Model* model,ConstDataHandle model_handle); - void CreateLoadsBedSlopeCompute(DataSet** ploads, Model* model, ConstDataHandle model_handle); - void CreateParametersBedSlopeCompute(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount); + /*slope compute*/ + void CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle); + void CreateConstraintsSlopeCompute(DataSet** pconstraints,Model* model,ConstDataHandle model_handle); + void CreateLoadsSlopeCompute(DataSet** ploads, Model* model, ConstDataHandle model_handle); + void CreateParametersSlopeCompute(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount); /*control:*/ Index: /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp (revision 357) +++ /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateConstraintsSlopeCompute.cpp (revision 358) @@ -1,8 +1,8 @@ /* - * CreateConstraintsBedSlopeCompute.c: + * CreateConstraintsSlopeCompute.c: */ #undef __FUNCT__ -#define __FUNCT__ "CreateConstraintsBedSlopeCompute" +#define __FUNCT__ "CreateConstraintsSlopeCompute" #include "../../DataSet/DataSet.h" @@ -14,5 +14,5 @@ -void CreateConstraintsBedSlopeCompute(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){ +void CreateConstraintsSlopeCompute(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){ Index: /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 357) +++ /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 358) @@ -1,8 +1,8 @@ /* - * CreateElementsNodesAndMaterialsBedSlopeCompute.c: + * CreateElementsNodesAndMaterialsSlopeCompute.c: */ #undef __FUNCT__ -#define __FUNCT__ "CreateElementsNodesAndMaterialsBedSlopeCompute" +#define __FUNCT__ "CreateElementsNodesAndMaterialsSlopeCompute" #include "../../DataSet/DataSet.h" @@ -15,6 +15,13 @@ -void CreateElementsNodesAndMaterialsBedSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){ - +void CreateElementsNodesAndMaterialsSlopeCompute(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: */ @@ -23,4 +30,94 @@ DataSet* materials = NULL; + /*Objects: */ + Node* node = NULL; + Tria* tria = NULL; + Penta* penta = NULL; + + int 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 + + /*tria constructor input: */ + int tria_id; + int tria_mid; + int tria_mparid; + int tria_g[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;/*!scaling ratio for velocities*/ + double tria_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ + double tria_viscosity_overshoot; + + /*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; + + /* 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; + 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()); @@ -31,4 +128,351 @@ if (!model->isstokes)goto cleanup_and_return; + /*Get analysis_type: */ + analysis_type=AnalysisTypeAsEnum(model->analysis_type); + + /*Width of elements: */ + if(strcmp(model->meshtype,"2d")==0){ + elements_width=3; //tria elements + } + else{ + elements_width=6; //penta elements + } + + #ifdef _PARALLEL_ + /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ + if(strcmp(model->meshtype,"2d")==0){ + /*load elements: */ + ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); + } + else{ + /*load elements2d: */ + 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->elements); + 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: */ + + /*2d mesh: */ + if (strcmp(model->meshtype,"2d")==0){ + + /*Fetch data needed: */ + ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); + ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat"); + ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat"); + + for (i=0;inumberofelements;i++){ + + #ifdef _PARALLEL_ + /*!All elements have been partitioned above, only create elements for this CPU: */ + if(my_rank==epart[i]){ + #endif + + + /*ids: */ + tria_id=i+1; //matlab indexing. + tria_mid=-1; //no materials + tria_mparid=-1; //no materials + + /*vertices offsets: */ + tria_g[0]=(int)*(model->elements+elements_width*i+0); + tria_g[1]=(int)*(model->elements+elements_width*i+1); + tria_g[2]=(int)*(model->elements+elements_width*i+2); + + /*surface and bed:*/ + tria_s[0]=*(model->surface+ ((int)*(model->elements+elements_width*i+0)-1)); + tria_s[1]=*(model->surface+ ((int)*(model->elements+elements_width*i+1)-1)); + tria_s[2]=*(model->surface+ ((int)*(model->elements+elements_width*i+2)-1)); + + tria_b[0]=*(model->bed+ ((int)*(model->elements+elements_width*i+0)-1)); + tria_b[1]=*(model->bed+ ((int)*(model->elements+elements_width*i+1)-1)); + tria_b[2]=*(model->bed+ ((int)*(model->elements+elements_width*i+2)-1)); + + /*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); + + /*Add tria element to elements dataset: */ + elements->AddObject(tria); + + #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; + #endif + + #ifdef _PARALLEL_ + }//if(my_rank==epart[i]) + #endif + + }//for (i=0;ielements); + xfree((void**)&model->surface); + xfree((void**)&model->bed); + + } + else{ // if (strcmp(meshtype,"2d")==0) + + /*Fetch data needed: */ + ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","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->elementonbed,NULL,NULL,model_handle,"elementonbed","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=-1; //no materials + penta_mparid=-1; //no materials + + + /*vertices,thickness,surface,bed and drag: */ + for(j=0;j<6;j++){ + penta_g[j]=(int)*(model->elements+elements_width*i+j); + 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)); + } + + /*diverse: */ + penta_onbed=(int)*(model->elementonbed+i); + + /*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); + + /*Add penta element to elements dataset: */ + elements->AddObject(penta); + + + #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->surface); + xfree((void**)&model->bed); + xfree((void**)&model->elementonbed); + + } //if (strcmp(meshtype,"2d")==0) + + #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); + + 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: @@ -37,3 +481,15 @@ *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/SlopeCompute/CreateLoadsSlopeCompute.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp (revision 357) +++ /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateLoadsSlopeCompute.cpp (revision 358) @@ -1,7 +1,7 @@ -/*! \file CreateLoadsBedSlopeCompute.c: +/*! \file CreateLoadsSlopeCompute.c: */ #undef __FUNCT__ -#define __FUNCT__ "CreateLoadsBedSlopeCompute" +#define __FUNCT__ "CreateLoadsSlopeCompute" #include "../../DataSet/DataSet.h" @@ -14,5 +14,5 @@ -void CreateLoadsBedSlopeCompute(DataSet** ploads, Model* model,ConstDataHandle model_handle){ +void CreateLoadsSlopeCompute(DataSet** ploads, Model* model,ConstDataHandle model_handle){ DataSet* loads = NULL; Index: /issm/trunk/src/c/SlopeExtrudex/SlopeExtrudex.cpp =================================================================== --- /issm/trunk/src/c/SlopeExtrudex/SlopeExtrudex.cpp (revision 358) +++ /issm/trunk/src/c/SlopeExtrudex/SlopeExtrudex.cpp (revision 358) @@ -0,0 +1,36 @@ +/*!\file SlopeExtrudex + * \brief: vertical velocity extrusion + */ + +#include "./SlopeExtrudex.h" + +#undef __FUNCT__ +#define __FUNCT__ "SlopeExtrudex" + +#include "../shared/shared.h" +#include "../include/macros.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void SlopeExtrudex( Vec sg, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials){ + + double* sg_serial=NULL; + + /*First, get elements and nodes configured: */ + elements->Configure(elements,loads, nodes, materials); + nodes->Configure(elements,loads, nodes, materials); + + /*Serialize velcoity: */ + VecToMPISerial(&sg_serial,sg); + + /*Extrude velocity vertically: */ + elements->SlopeExtrude(sg,sg_serial); + + /*Assemble vector: */ + VecAssemblyBegin(sg); + VecAssemblyEnd(sg); + + /*Free ressources:*/ + xfree((void**)&sg_serial); + +} Index: /issm/trunk/src/c/SlopeExtrudex/SlopeExtrudex.h =================================================================== --- /issm/trunk/src/c/SlopeExtrudex/SlopeExtrudex.h (revision 358) +++ /issm/trunk/src/c/SlopeExtrudex/SlopeExtrudex.h (revision 358) @@ -0,0 +1,14 @@ +/*!\file: SlopeExtrudex.h + * \brief header file for velocity extrusion + */ + +#ifndef _SLOPEEXTRUDEX_H +#define _SLOPEEXTRUDEX_H + +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void SlopeExtrudex( Vec sg, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials); + +#endif /* _SLOPEEXTRUDEX_H */ + Index: /issm/trunk/src/c/issm.h =================================================================== --- /issm/trunk/src/c/issm.h (revision 357) +++ /issm/trunk/src/c/issm.h (revision 358) @@ -47,4 +47,5 @@ #include "./GriddataMeshToGridx/GriddataMeshToGridx.h" #include "./ComputePressurex/ComputePressurex.h" +#include "./SlopeExtrudex/SlopeExtrudex.h" Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 357) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 358) @@ -292,4 +292,13 @@ } + else if ( + (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || + (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) || + (analysis_type==BedSlopeComputeXAnalysisEnum()) || + (analysis_type==BedSlopeComputeYAnalysisEnum()) + ){ + + CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type); + } else{ throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet")); @@ -703,4 +712,13 @@ CreatePVectorDiagnosticVert( pg,inputs,analysis_type); } + else if ( + (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || + (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) || + (analysis_type==BedSlopeComputeXAnalysisEnum()) || + (analysis_type==BedSlopeComputeYAnalysisEnum()) + ){ + + CreatePVectorSlopeCompute( pg,inputs,analysis_type); + } else{ throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet")); @@ -1575,4 +1593,55 @@ #undef __FUNCT__ +#define __FUNCT__ "Penta::SlopeExtrude" +void Penta::SlopeExtrude(Vec sg,double* sg_serial){ + + /* node data: */ + const int numgrids=6; + const int NDOF1=1; + const int numdofs=NDOF1*numgrids; + int doflist[numdofs]; + int nodedof; + int numberofdofspernode; + + Node* node=NULL; + int i; + double slope; + + + if(onbed==1){ + + GetDofList(&doflist[0],&numberofdofspernode); + + if(numberofdofspernode!=1)throw ErrorException(__FUNCT__," slope can only be extruded on 1 dof per node"); + + /*For each node on the base of this penta, we grab the slope. Once we know the slope, we follow the upper nodes, + * inserting the same slope value into sg, until we reach the surface: */ + + for(i=0;i<3;i++){ + + node=nodes[i]; //base nodes + + /*get velocity for this base node: */ + slope=sg_serial[doflist[i]]; + + //go throsgn all nodes which sit on top of this node, until we reach the surface, + //and plsg slope in sg + for(;;){ + + node->GetDofList(&nodedof,&numberofdofspernode); + VecSetValues(sg,1,&nodedof,&slope,INSERT_VALUES); + + if (node->IsOnSurface())break; + /*get next node: */ + node=node->GetUpperNode(); + } + } + + } + +} + + +#undef __FUNCT__ #define __FUNCT__ "Penta:GetB_vert" void Penta::GetB_vert(double* B, double* xyz_list, double* gauss_l1l2l3l4){ @@ -1770,5 +1839,5 @@ double pressure[numgrids]; double rho_ice,g; - double xyz_list[numgrids][3]; + double xyz_list[numgrids][3]; /*Get node data: */ @@ -1792,2 +1861,41 @@ } + + +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreateKMatrixSlopeCompute" + +void Penta::CreateKMatrixSlopeCompute(Mat Kgg,void* inputs,int analysis_type){ + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreateKMatrix(Kgg,inputs, analysis_type); + delete tria; + return; + +} + +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreatePVectorSlopeCompute" + +void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type){ + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreatePVector(pg,inputs, analysis_type); + delete tria; + return; +} + Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 357) +++ /issm/trunk/src/c/objects/Penta.h (revision 358) @@ -110,5 +110,9 @@ void GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4); void VelocityExtrude(Vec ug,double* ug_serial); + void SlopeExtrude(Vec sg,double* sg_serial); void ComputePressure(Vec p_g); + void CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type); + void CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type); + }; Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 357) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 358) @@ -274,4 +274,13 @@ } + else if ( + (analysis_type==SurfaceSlopeComputeXAnalysisEnum()) || + (analysis_type==SurfaceSlopeComputeYAnalysisEnum()) || + (analysis_type==BedSlopeComputeXAnalysisEnum()) || + (analysis_type==BedSlopeComputeYAnalysisEnum()) + ){ + + CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type); + } else{ throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet")); @@ -666,4 +675,86 @@ #undef __FUNCT__ +#define __FUNCT__ "Tria::CreateKMatrixSlopeCompute" + +void Tria::CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type){ + + /* local declarations */ + int i,j; + + /* node data: */ + const int numgrids=3; + const int NDOF1=1; + const int numdofs=NDOF1*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdofs]; + int numberofdofspernode; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + double gauss_weight; + double gauss_l1l2l3[3]; + + /* matrices: */ + double L[1][3]; + double DL_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 Jdet; + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Set Ke_gg to 0: */ + for(i=0;iFindParam((void*)&u_g_initial,"u_g"); femmodel.parameters->FindParam((void*)&numberofnodes,"numberofnodes"); - + inputs=new ParameterInputs; inputs->Add("velocity",u_g_initial,3,numberofnodes); Index: /issm/trunk/src/c/shared/Elements/ResolvePointers.cpp =================================================================== --- /issm/trunk/src/c/shared/Elements/ResolvePointers.cpp (revision 357) +++ /issm/trunk/src/c/shared/Elements/ResolvePointers.cpp (revision 358) @@ -29,4 +29,6 @@ for(i=0;i