Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 3416) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 3417) @@ -7,6 +7,6 @@ -#ifndef DATASET_H_ -#define DATASET_H_ +#ifndef _DATASET_H_ +#define _DATASET_H_ #include Index: /issm/trunk/src/c/Dofx/Dofx.cpp =================================================================== --- /issm/trunk/src/c/Dofx/Dofx.cpp (revision 3416) +++ /issm/trunk/src/c/Dofx/Dofx.cpp (revision 3417) @@ -10,5 +10,5 @@ #include "../EnumDefinitions/EnumDefinitions.h" -int Dofx( DofVec** ppartition, DofVec** ptpartition,DataSet* elements,DataSet* nodes, DataSet* params) { +int Dofx( DofVec** ppartition, DofVec** ptpartition, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* params) { int noerr=1; @@ -21,4 +21,5 @@ /*intermediary: */ int numberofnodes; + int numberofvertices; int numberofdofspernode; @@ -31,5 +32,8 @@ tpartition=new DofVec("tpartition"); - /*First, recover number of grids from parameters: */ + /*First, recover number of vertices and nodes from parameters: */ + found=params->FindParam(&numberofvertices,"numberofvertices"); + if(!found)ISSMERROR(" could not find numberofvertices in parameters"); + found=params->FindParam(&numberofnodes,"numberofnodes"); if(!found)ISSMERROR(" could not find numberofnodes in parameters"); @@ -41,4 +45,5 @@ /*Ensure that only for each cpu, the partition border nodes only will be taken into account once * across the cluster. To do so, we flag all the clone nodes: */ + vertices->FlagClones(numberofvertices); nodes->FlagClones(numberofnodes); Index: /issm/trunk/src/c/Dofx/Dofx.h =================================================================== --- /issm/trunk/src/c/Dofx/Dofx.h (revision 3416) +++ /issm/trunk/src/c/Dofx/Dofx.h (revision 3417) @@ -9,5 +9,5 @@ /* local prototypes: */ -int Dofx( DofVec** partition, DofVec** ptpartition,DataSet* elements,DataSet* nodesin,DataSet* params); +int Dofx( DofVec** partition, DofVec** ptpartition,DataSet* elements,DataSet* nodesin, DataSet* verticesin, DataSet* params); #endif /* _DOFX_H */ Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 3416) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 3417) @@ -15,4 +15,5 @@ int ParametersEnum(void){ return 106; } int ResultsEnum(void){ return 107; } +int VerticesEnum(void){ return 108; } /*ANALYSIS TYPES: */ @@ -71,10 +72,14 @@ int ElementEnum(void){ return 410; } int TriaEnum(void){ return 411; } -int ElementPropertiesEnum(void){ return 412; } -int PentaEnum(void){ return 413; } -int SingEnum(void){ return 414; } -int BeamEnum(void){ return 415; } +int ElementPropertiesEnum(void){ return 412; } +int NodePropertiesEnum(void){ return 413; } +int PentaEnum(void){ return 414; } +int SingEnum(void){ return 415; } +int BeamEnum(void){ return 416; } +int DofIndexingEnum(void){ return 417; } + /*Grids: */ -int NodeEnum(void){ return 420; } +int VertexEnum(void){ return 420; } +int NodeEnum(void){ return 421; } /*Loads: */ int LoadEnum(void){ return 430; } Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 3416) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 3417) @@ -16,4 +16,5 @@ int ParametersEnum(void); int ResultsEnum(void); +int VerticesEnum(void); /*ANALYSIS TYPES: */ @@ -73,9 +74,12 @@ int TriaEnum(void); int ElementPropertiesEnum(void); +int NodePropertiesEnum(void); int PentaEnum(void); int SingEnum(void); int BeamEnum(void); +int DofIndexingEnum(void); /*Grids: */ int NodeEnum(void); +int VertexEnum(void); /*Loads: */ int LoadEnum(void); Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 3416) +++ /issm/trunk/src/c/Makefile.am (revision 3417) @@ -52,10 +52,14 @@ ./objects/DakotaPlugin.h\ ./objects/DakotaPlugin.cpp\ - ./objects/Node.h\ - ./objects/Node.cpp\ ./objects/Vertex.h\ ./objects/Vertex.cpp\ ./objects/Hook.h\ ./objects/Hook.cpp\ + ./objects/DofIndexing.h\ + ./objects/DofIndexing.cpp\ + ./objects/NodeProperties.h\ + ./objects/NodeProperties.cpp\ + ./objects/Node.h\ + ./objects/Node.cpp\ ./objects/Result.h\ ./objects/Result.cpp\ @@ -201,4 +205,5 @@ ./ModelProcessorx/IoModel.h\ ./ModelProcessorx/IoModel.cpp\ + ./ModelProcessorx/Partitioning.cpp\ ./ModelProcessorx/CreateDataSets.cpp\ ./ModelProcessorx/CreateParameters.cpp\ @@ -446,4 +451,8 @@ ./objects/DakotaPlugin.h\ ./objects/DakotaPlugin.cpp\ + ./objects/DofIndexing.h\ + ./objects/DofIndexing.cpp\ + ./objects/NodeProperties.h\ + ./objects/NodeProperties.cpp\ ./objects/Node.h\ ./objects/Node.cpp\ @@ -592,4 +601,5 @@ ./ModelProcessorx/IoModel.h\ ./ModelProcessorx/IoModel.cpp\ + ./ModelProcessorx/Partitioning.cpp\ ./ModelProcessorx/CreateDataSets.cpp\ ./ModelProcessorx/CreateParameters.cpp\ Index: /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp (revision 3416) +++ /issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp (revision 3417) @@ -13,5 +13,5 @@ -void CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ +void CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ @@ -26,8 +26,11 @@ DataSet* elements = NULL; DataSet* nodes = NULL; + DataSet* vertices = NULL; DataSet* materials = NULL; /*Objects: */ Node* node = NULL; + NodeProperties node_properties; + Vertex* vertex = NULL; Tria* tria = NULL; Penta* penta = NULL; @@ -96,7 +99,6 @@ /* node constructor input: */ int node_id; - int node_partitionborder=0; - double node_x[3]; - double node_sigma; + int vertex_id; + int partitionborder=0; int node_onbed; int node_onsurface; @@ -121,4 +123,5 @@ elements = new DataSet(ElementsEnum()); nodes = new DataSet(NodesEnum()); + vertices = new DataSet(VerticesEnum()); materials = new DataSet(MaterialsEnum()); @@ -436,28 +439,29 @@ #endif + #ifdef _PARALLEL_ + if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border + partitionborder=1; + } + else{ + partitionborder=0; + } + #else + partitionborder=0; + #endif + + /*create vertex: */ + vertex_id=i+1; + vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]),partitionborder); + vertices->AddObject(vertex); + + /*create node: */ 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]=iomodel->x[i]; - node_x[1]=iomodel->y[i]; - node_x[2]=iomodel->z[i]; - node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]); - - node_onbed=(int)iomodel->gridonbed[i]; - node_onsurface=(int)iomodel->gridonsurface[i]; - node_onshelf=(int)iomodel->gridoniceshelf[i]; - node_onsheet=(int)iomodel->gridonicesheet[i]; + node_properties.SetProperties( + (int)iomodel->gridonbed[i], + (int)iomodel->gridonsurface[i], + (int)iomodel->gridoniceshelf[i], + (int)iomodel->gridonicesheet[i]); if (strcmp(iomodel->meshtype,"3d")==0){ @@ -475,5 +479,6 @@ /*Create node using its constructor: */ - node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet); + node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties); + /*set single point constraints.: */ @@ -497,4 +502,5 @@ * datasets, it will not be redone: */ elements->Presort(); + vertices->Presort(); nodes->Presort(); materials->Presort(); @@ -531,4 +537,5 @@ *pelements=elements; *pnodes=nodes; + *pvertices=vertices; *pmaterials=materials; Index: /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp (revision 3416) +++ /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp (revision 3417) @@ -13,5 +13,5 @@ -void CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ +void CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ @@ -26,4 +26,5 @@ DataSet* elements = NULL; DataSet* nodes = NULL; + DataSet* vertices = NULL; DataSet* materials = NULL; ElementProperties* tria_properties=NULL; @@ -32,4 +33,5 @@ /*Objects: */ Node* node = NULL; + Vertex* vertex = NULL; Tria* tria = NULL; Penta* penta = NULL; @@ -96,7 +98,6 @@ /* node constructor input: */ int node_id; - int node_partitionborder=0; - double node_x[3]; - double node_sigma; + int vertex_id; + int partitionborder=0; int node_onbed; int node_onsurface; @@ -121,4 +122,5 @@ elements = new DataSet(ElementsEnum()); nodes = new DataSet(NodesEnum()); + vertices = new DataSet(VerticesEnum()); materials = new DataSet(MaterialsEnum()); @@ -435,8 +437,12 @@ #endif + /*create vertex: */ + vertex_id=i+1; + vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i])); + + vertices->AddObject(vertex); + node_id=i+1; //matlab indexing - - #ifdef _PARALLEL_ if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border @@ -450,13 +456,10 @@ #endif - node_x[0]=iomodel->x[i]; - node_x[1]=iomodel->y[i]; - node_x[2]=iomodel->z[i]; - node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]); - - node_onbed=(int)iomodel->gridonbed[i]; - node_onsurface=(int)iomodel->gridonsurface[i]; - node_onshelf=(int)iomodel->gridoniceshelf[i]; - node_onsheet=(int)iomodel->gridonicesheet[i]; + node_properties.SetProperties( + (int)iomodel->gridonbed[i], + (int)iomodel->gridonsurface[i], + (int)iomodel->gridoniceshelf[i], + (int)iomodel->gridonicesheet[i]); + if (strcmp(iomodel->meshtype,"3d")==0){ @@ -474,5 +477,5 @@ /*Create node using its constructor: */ - node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet); + node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties); /*set single point constraints.: */ @@ -497,4 +500,5 @@ elements->Presort(); nodes->Presort(); + vertices->Presort(); materials->Presort(); @@ -530,4 +534,5 @@ *pelements=elements; *pnodes=nodes; + *pvertices=vertices; *pmaterials=materials; Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 3416) +++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 3417) @@ -15,5 +15,5 @@ -void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){ +void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){ /*create parameters common to all solutions: */ @@ -26,5 +26,5 @@ if (iomodel->sub_analysis_type==HorizAnalysisEnum()){ - CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes,pmaterials, iomodel,iomodel_handle); + CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle); CreateConstraintsDiagnosticHoriz(pconstraints,iomodel,iomodel_handle); CreateLoadsDiagnosticHoriz(ploads,iomodel,iomodel_handle); @@ -34,5 +34,5 @@ else if (iomodel->sub_analysis_type==VertAnalysisEnum()){ - CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes,pmaterials, iomodel,iomodel_handle); + CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle); CreateConstraintsDiagnosticVert(pconstraints,iomodel,iomodel_handle); CreateLoadsDiagnosticVert(ploads,iomodel,iomodel_handle); @@ -41,5 +41,5 @@ else if (iomodel->sub_analysis_type==StokesAnalysisEnum()){ - CreateElementsNodesAndMaterialsDiagnosticStokes(pelements,pnodes,pmaterials, iomodel,iomodel_handle); + CreateElementsNodesAndMaterialsDiagnosticStokes(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle); CreateConstraintsDiagnosticStokes(pconstraints,iomodel,iomodel_handle); CreateLoadsDiagnosticStokes(ploads,iomodel,iomodel_handle); @@ -48,5 +48,5 @@ else if (iomodel->sub_analysis_type==HutterAnalysisEnum()){ - CreateElementsNodesAndMaterialsDiagnosticHutter(pelements,pnodes,pmaterials, iomodel,iomodel_handle); + CreateElementsNodesAndMaterialsDiagnosticHutter(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); CreateConstraintsDiagnosticHutter(pconstraints,iomodel,iomodel_handle); CreateLoadsDiagnosticHutter(ploads,iomodel,iomodel_handle); @@ -56,5 +56,5 @@ else if (iomodel->analysis_type==SlopecomputeAnalysisEnum()){ - CreateElementsNodesAndMaterialsSlopeCompute(pelements,pnodes,pmaterials, iomodel,iomodel_handle); + CreateElementsNodesAndMaterialsSlopeCompute(pelements,pnodes, pvertices,pmaterials, iomodel,iomodel_handle); CreateConstraintsSlopeCompute(pconstraints,iomodel,iomodel_handle); CreateLoadsSlopeCompute(ploads,iomodel,iomodel_handle); @@ -63,5 +63,5 @@ else if (iomodel->analysis_type==ThermalAnalysisEnum()){ - CreateElementsNodesAndMaterialsThermal(pelements,pnodes,pmaterials, iomodel,iomodel_handle); + CreateElementsNodesAndMaterialsThermal(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); CreateConstraintsThermal(pconstraints,iomodel,iomodel_handle); CreateLoadsThermal(ploads,iomodel,iomodel_handle); @@ -71,5 +71,5 @@ else if (iomodel->analysis_type==MeltingAnalysisEnum()){ - CreateElementsNodesAndMaterialsMelting(pelements,pnodes,pmaterials, iomodel,iomodel_handle); + CreateElementsNodesAndMaterialsMelting(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); CreateConstraintsMelting(pconstraints,iomodel,iomodel_handle); CreateLoadsMelting(ploads,iomodel,iomodel_handle); @@ -78,5 +78,5 @@ else if (iomodel->analysis_type==PrognosticAnalysisEnum()){ - CreateElementsNodesAndMaterialsPrognostic(pelements,pnodes,pmaterials, iomodel,iomodel_handle); + CreateElementsNodesAndMaterialsPrognostic(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); CreateConstraintsPrognostic(pconstraints,iomodel,iomodel_handle); CreateLoadsPrognostic(ploads,iomodel,iomodel_handle); @@ -86,5 +86,5 @@ else if (iomodel->analysis_type==Prognostic2AnalysisEnum()){ - CreateElementsNodesAndMaterialsPrognostic2(pelements,pnodes,pmaterials, iomodel,iomodel_handle); + CreateElementsNodesAndMaterialsPrognostic2(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); CreateConstraintsPrognostic2(pconstraints,iomodel,iomodel_handle); CreateLoadsPrognostic2(ploads,iomodel,iomodel_handle); @@ -94,5 +94,5 @@ else if (iomodel->analysis_type==BalancedthicknessAnalysisEnum()){ - CreateElementsNodesAndMaterialsBalancedthickness(pelements,pnodes,pmaterials, iomodel,iomodel_handle); + CreateElementsNodesAndMaterialsBalancedthickness(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); CreateConstraintsBalancedthickness(pconstraints,iomodel,iomodel_handle); CreateLoadsBalancedthickness(ploads,iomodel,iomodel_handle); @@ -102,5 +102,5 @@ else if (iomodel->analysis_type==BalancedvelocitiesAnalysisEnum()){ - CreateElementsNodesAndMaterialsBalancedvelocities(pelements,pnodes,pmaterials, iomodel,iomodel_handle); + CreateElementsNodesAndMaterialsBalancedvelocities(pelements,pnodes,pvertices,pmaterials, iomodel,iomodel_handle); CreateConstraintsBalancedvelocities(pconstraints,iomodel,iomodel_handle); CreateLoadsBalancedvelocities(ploads,iomodel,iomodel_handle); Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 3416) +++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 3417) @@ -238,9 +238,15 @@ parameters->AddObject(param); - /*numberofgrids: */ + /*numberofvertices: */ + count++; + param= new Param(count,"numberofvertices",DOUBLE); + param->SetDouble(iomodel->numberofvertices); + parameters->AddObject(param); + + /*numberofnodes: */ count++; param= new Param(count,"numberofnodes",DOUBLE); if (iomodel->analysis_type==Prognostic2AnalysisEnum()) param->SetDouble(3*iomodel->numberofelements); - else param->SetDouble(iomodel->numberofnodes); + else param->SetDouble(iomodel->numberofvertices); parameters->AddObject(param); Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 3416) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 3417) @@ -13,178 +13,34 @@ -void CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ - - - /*output: int* epart, int* my_grids, double* my_bordergrids*/ +void CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ int i,j,k,n; - extern int my_rank; - extern int num_procs; /*DataSets: */ DataSet* elements = NULL; DataSet* nodes = NULL; + DataSet* vertices = NULL; DataSet* materials = NULL; /*Objects: */ Node* node = NULL; + Vertex* vertex = NULL; Tria* tria = NULL; Penta* penta = NULL; Matice* matice = NULL; Matpar* matpar = NULL; - ElementProperties* tria_properties=NULL; - ElementProperties* penta_properties=NULL; - - /*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; - - /*tria constructor input: */ - int tria_id; - int tria_matice_id; - int tria_matpar_id; - int tria_numpar_id; - 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; - bool tria_onwater; - - /*matice constructor input: */ - int matice_mid; - double matice_B; - double matice_n; - - /*penta constructor input: */ - int penta_id; - int penta_matice_id; - int penta_matpar_id; - int penta_numpar_id; - int penta_node_ids[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; - int penta_collapse; - double penta_melting[6]; - double penta_accumulation[6]; - bool penta_onwater; - - /*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]; - double node_sigma; - int node_onbed; - int node_onsurface; - int node_onsheet; - int node_onshelf; - 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 el1,el2; - + + /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ + if (!iomodel->ismacayealpattyn)goto cleanup_and_return; + /*First create the elements, nodes and material properties: */ elements = new DataSet(ElementsEnum()); nodes = new DataSet(NodesEnum()); + vertices = new DataSet(VerticesEnum()); materials = new DataSet(MaterialsEnum()); - - /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ - if (!iomodel->ismacayealpattyn)goto cleanup_and_return; - - /*Width of elements: */ - if(strcmp(iomodel->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(iomodel->meshtype,"2d")==0){ - /*load elements: */ - IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); - } - else{ - /*load elements2d: */ - IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d"); - } - - MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs); - - /*Free elements and elements2d: */ - xfree((void**)&iomodel->elements); - xfree((void**)&iomodel->elements2d); - - - /*Deal with rifts, they have to be included into one partition only, not several: */ - IoModelFetchData(&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo"); - - for(i=0;inumrifts;i++){ - el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2)-1; //matlab indexing to c indexing - el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+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**)&iomodel->riftinfo); - - /*Used later on: */ - my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int)); - #endif - - - - /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ + /*Partition elements and vertices and nodes: */ + Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle); /*2d mesh: */ @@ -209,107 +65,16 @@ 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 - - if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements. + if(iomodel->my_elements[i]){ - /*ids: */ - tria_id=i+1; //matlab indexing. - tria_matice_id=i+1; //refers to the corresponding material property card - tria_matpar_id=iomodel->numberofelements+1;//refers to the corresponding parmat property card - tria_numpar_id=1; - - /*vertices ids: */ - tria_node_ids[0]=(int)*(iomodel->elements+elements_width*i+0); - tria_node_ids[1]=(int)*(iomodel->elements+elements_width*i+1); - tria_node_ids[2]=(int)*(iomodel->elements+elements_width*i+2); - - /*thickness,surface and bed:*/ - tria_h[0]= *(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing. - tria_h[1]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+1)-1)); - tria_h[2]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+2)-1)) ; - - tria_s[0]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+0)-1)); - tria_s[1]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+1)-1)); - tria_s[2]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+2)-1)); - - tria_b[0]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+0)-1)); - tria_b[1]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+1)-1)); - tria_b[2]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+2)-1)); - - /*basal drag:*/ - tria_friction_type=(int)iomodel->drag_type; - - tria_k[0]=*(iomodel->drag+ ((int)*(iomodel->elements+elements_width*i+0)-1)); - tria_k[1]=*(iomodel->drag+ ((int)*(iomodel->elements+elements_width*i+1)-1)); - tria_k[2]=*(iomodel->drag+ ((int)*(iomodel->elements+elements_width*i+2)-1)); - - tria_p=iomodel->p[i]; - tria_q=iomodel->q[i]; - - /*meling and accumulation*/ - tria_melting[0]=*(iomodel->melting+ ((int)*(iomodel->elements+elements_width*i+0)-1)); - tria_melting[1]=*(iomodel->melting+ ((int)*(iomodel->elements+elements_width*i+1)-1)); - tria_melting[2]=*(iomodel->melting+ ((int)*(iomodel->elements+elements_width*i+2)-1)); - - tria_accumulation[0]=*(iomodel->accumulation+ ((int)*(iomodel->elements+elements_width*i+0)-1)); - tria_accumulation[1]=*(iomodel->accumulation+ ((int)*(iomodel->elements+elements_width*i+1)-1)); - tria_accumulation[2]=*(iomodel->accumulation+ ((int)*(iomodel->elements+elements_width*i+2)-1)); - - /*element on iceshelf, water?:*/ - tria_shelf=(int)*(iomodel->elementoniceshelf+i); - tria_onwater=(bool)*(iomodel->elementonwater+i); - - /*Create properties: */ - tria_properties=new ElementProperties(3,tria_h, tria_s, tria_b, tria_k, tria_melting, tria_accumulation, NULL, - tria_friction_type, tria_p, tria_q, tria_shelf, UNDEF,tria_onwater, UNDEF,UNDEF,UNDEF); - - /*Create tria element using its constructor:*/ - tria=new Tria(tria_id, tria_node_ids, tria_matice_id, tria_matpar_id, tria_numpar_id, tria_properties); - - /*delete properties: */ - delete tria_properties; - - /*Add tria element to elements dataset: */ - elements->AddObject(tria); - - /*Deal with material property card: */ - matice_mid=i+1; //same as the material id from the geom2 elements. - - /*Average B over 3 grid elements: */ - B_avg=0; - for(j=0;j<3;j++){ - B_avg+=*(iomodel->B+((int)*(iomodel->elements+elements_width*i+j)-1)); + if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements. + + /*Create and add tria element to elements dataset: */ + elements->AddObject(new Tria(i,iomodel)); + + /*Create and add material property to materials dataset: */ + materials->AddObject(new Matice(i,iomodel,3)); } - B_avg=B_avg/3; - matice_B=B_avg; - matice_n=(double)*(iomodel->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); - - } //if(!MacAyealFormulationEnum) - - #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)*(iomodel->elements+elements_width*i+0)-1]=1; - my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1; - my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1; - #endif - - #ifdef _PARALLEL_ - }//if(my_rank==epart[i]) - #endif - + } }//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 - - if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum() | *(iomodel->elements_type+2*i+0)==PattynFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements. - - /*name and id: */ - penta_id=i+1; //matlab indexing. - penta_matice_id=i+1; //refers to the corresponding material property card - penta_matpar_id=iomodel->numberofelements+1;//refers to the corresponding parmat property card - penta_numpar_id=1; - - /*vertices,thickness,surface,bed and drag: */ - for(j=0;j<6;j++){ - penta_node_ids[j]=(int)*(iomodel->elements+elements_width*i+j); - penta_h[j]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+j)-1)); - penta_s[j]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+j)-1)); - penta_b[j]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+j)-1)); - penta_k[j]=*(iomodel->drag+ ((int)*(iomodel->elements+elements_width*i+j)-1)); - penta_melting[j]=*(iomodel->melting+ ((int)*(iomodel->elements+elements_width*i+j)-1)); - penta_accumulation[j]=*(iomodel->accumulation+ ((int)*(iomodel->elements+elements_width*i+j)-1)); + if(my_elements[i]){ + if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum() | *(iomodel->elements_type+2*i+0)==PattynFormulationEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements. + /*Create and add penta element to elements dataset: */ + elements->AddObject(new Penta(i,iomodel)); + + /*Create and add material property to materials dataset: */ + materials->AddObject(new Matice(i,iomodel,6)); + } - - /*basal drag:*/ - penta_friction_type=(int)iomodel->drag_type; - - penta_p=iomodel->p[i]; - penta_q=iomodel->q[i]; - - /*diverse: */ - penta_shelf=(int)*(iomodel->elementoniceshelf+i); - penta_onbed=(int)*(iomodel->elementonbed+i); - penta_onsurface=(int)*(iomodel->elementonsurface+i); - penta_onwater=(bool)*(iomodel->elementonwater+i); - - if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base. - penta_collapse=1; - } - else{ - penta_collapse=0; - } - - /*Create element properties: */ - penta_properties=new ElementProperties(6,penta_h, penta_s, penta_b, penta_k, penta_melting, penta_accumulation, NULL, penta_friction_type, penta_p, penta_q, penta_shelf, penta_onbed, penta_onwater, penta_onsurface, penta_collapse, UNDEF); - - /*Create Penta using its constructor:*/ - penta= new Penta(penta_id,penta_node_ids, penta_matice_id, penta_matpar_id, penta_numpar_id, penta_properties); - - /*delete properties: */ - delete penta_properties; - - /*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+=*(iomodel->B+((int)*(iomodel->elements+elements_width*i+j)-1)); - } - B_avg=B_avg/6; - matice_B= B_avg; - matice_n=(double)*(iomodel->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)*(iomodel->elements+elements_width*i+0)-1]=1; - my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1; - my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1; - my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1; - my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1; - my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1; - #endif - - #ifdef _PARALLEL_ - }//if(my_rank==epart[i]) - #endif + }//if(my_elements[i]) }//for (i=0;inumberofelements+1; //put it at the end of the materials - matpar_g=iomodel->g; - matpar_rho_ice=iomodel->rho_ice; - matpar_rho_water=iomodel->rho_water; - matpar_thermalconductivity=iomodel->thermalconductivity; - matpar_heatcapacity=iomodel->heatcapacity; - matpar_latentheat=iomodel->latentheat; - matpar_beta=iomodel->beta; - matpar_meltingpoint=iomodel->meltingpoint; - matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; - matpar_thermal_exchange_velocity=iomodel->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(iomodel->numberofnodes); - - for (i=0;inumberofnodes;i++){ - if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); - } - VecAssemblyBegin(gridborder); - VecAssemblyEnd(gridborder); - - #ifdef _ISSM_DEBUG_ - VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD); - #endif - - VecToMPISerial(&my_bordergrids,gridborder); - - #ifdef _ISSM_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(iomodel->meshtype,"3d")==0){ - - /*Get penalties: */ - IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties"); - - if(iomodel->numpenalties){ - - iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int)); - #ifdef _SERIAL_ - for(i=0;inumpenalties;i++)iomodel->penaltypartitioning[i]=1; - #else - for(i=0;inumpenalties;i++)iomodel->penaltypartitioning[i]=-1; - - for(i=0;inumpenalties;i++){ - first_grid_index=(int)(*(iomodel->penalties+i*iomodel->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.:*/ - iomodel->penaltypartitioning[i]=1; - } - if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border - iomodel->penaltypartitioning[i]=0; - } - } - #endif - } - - /*Free penalties: */ - xfree((void**)&iomodel->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(iomodel->meshtype,"3d")==0){ - IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); - IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); - } + + /*Add new constrant material property tgo materials, at the end: */ + materials->AddObject(new Matpar(iomodel)); + + /*Create nodes and vertices: */ IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x"); IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y"); @@ -563,82 +162,23 @@ IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet"); IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); - - /*Get number of dofs per node: */ - DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->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 + if (strcmp(iomodel->meshtype,"3d")==0){ + IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); + IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); + } + + + for (i=0;inumberovertices;i++){ + + /*vertices and nodes (same number, as we are running continuous galerkin formulation: */ + if(iomodel->my_vertices[i]){ - - #ifdef _PARALLEL_ - if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border - node_partitionborder=1; + /*Add vertex to vertices dataset: */ + vertices->AddObject(new Vertex(i,iomodel)); + + /*Add node to nodes dataset: */ + nodes->AddObject(new Node(i,iomodel)); + } - else{ - node_partitionborder=0; - } - #else - node_partitionborder=0; - #endif - - node_x[0]=iomodel->x[i]; - node_x[1]=iomodel->y[i]; - node_x[2]=iomodel->z[i]; - node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]); - - node_onbed=(int)iomodel->gridonbed[i]; - node_onsurface=(int)iomodel->gridonsurface[i]; - node_onshelf=(int)iomodel->gridoniceshelf[i]; - node_onsheet=(int)iomodel->gridonicesheet[i]; - - if (strcmp(iomodel->meshtype,"3d")==0){ - if (isnan(iomodel->uppernodes[i])){ - node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves. - } - else{ - node_upper_node_id=(int)iomodel->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_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet); - - /*set single point constraints.: */ - if (strcmp(iomodel->meshtype,"3d")==0){ - /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */ - if (iomodel->deadgrids[i]){ - for(k=1;k<=node_numdofs;k++){ - node->FreezeDof(k); - } - } - } - if (iomodel->gridonhutter[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: */ @@ -655,23 +195,16 @@ xfree((void**)&iomodel->gridonicesheet); xfree((void**)&iomodel->gridoniceshelf); - - - /*Keep partitioning information into iomodel*/ - iomodel->epart=epart; - iomodel->my_grids=my_grids; - iomodel->my_bordergrids=my_bordergrids; - - /*Free ressources:*/ - #ifdef _PARALLEL_ - xfree((void**)&all_numgrids); - xfree((void**)&npart); - VecFree(&gridborder); - #endif - - cleanup_and_return: - + + /*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(); + vertices->Presort(); + materials->Presort(); + /*Assign output pointer: */ *pelements=elements; *pnodes=nodes; + *pvertices=vertices; *pmaterials=materials; Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (revision 3416) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (revision 3417) @@ -13,5 +13,5 @@ #include "../IoModel.h" -void CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ +void CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ @@ -26,8 +26,10 @@ DataSet* elements = NULL; DataSet* nodes = NULL; + DataSet* vertices = NULL; DataSet* materials = NULL; /*Objects: */ Node* node = NULL; + Vertex* vertex = NULL; Penta* penta = NULL; Matice* matice = NULL; @@ -88,7 +90,6 @@ /* node constructor input: */ int node_id; + int vertex_id; int node_partitionborder=0; - double node_x[3]; - double node_sigma; int node_onbed; int node_onsurface; @@ -113,4 +114,5 @@ elements = new DataSet(ElementsEnum()); nodes = new DataSet(NodesEnum()); + vertices = new DataSet(VerticesEnum()); materials = new DataSet(MaterialsEnum()); @@ -400,7 +402,11 @@ #endif + /*create vertex: */ + vertex_id=i+1; + vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i])); + + vertices->AddObject(vertex); + node_id=i+1; //matlab indexing - - #ifdef _PARALLEL_ @@ -415,14 +421,11 @@ #endif - node_x[0]=iomodel->x[i]; - node_x[1]=iomodel->y[i]; - node_x[2]=iomodel->z[i]; - node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]); - - node_onbed=(int)iomodel->gridonbed[i]; - node_onsurface=(int)iomodel->gridonsurface[i]; - node_onshelf=(int)iomodel->gridoniceshelf[i]; - node_onsheet=(int)iomodel->gridonicesheet[i]; - + node_properties.SetProperties( + (int)iomodel->gridonbed[i], + (int)iomodel->gridonsurface[i], + (int)iomodel->gridoniceshelf[i], + (int)iomodel->gridonicesheet[i]); + + if (strcmp(iomodel->meshtype,"3d")==0){ if (isnan(iomodel->uppernodes[i])){ @@ -439,5 +442,5 @@ /*Create node using its constructor: */ - node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet); + node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties); /*set single point constraints.: */ @@ -468,4 +471,5 @@ elements->Presort(); nodes->Presort(); + vertices->Presort(); materials->Presort(); @@ -503,4 +507,5 @@ *pelements=elements; *pnodes=nodes; + *pvertices=vertices; *pmaterials=materials; Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 3416) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 3417) @@ -14,5 +14,5 @@ -void CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ +void CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ @@ -27,8 +27,10 @@ DataSet* elements = NULL; DataSet* nodes = NULL; + DataSet* vertices = NULL; DataSet* materials = NULL; /*Objects: */ Node* node = NULL; + Vertex* vertex = NULL; Penta* penta = NULL; Matice* matice = NULL; @@ -88,7 +90,6 @@ /* node constructor input: */ int node_id; + int vertex_id; int node_partitionborder=0; - double node_x[3]; - double node_sigma; int node_onbed; int node_onsurface; @@ -123,4 +124,5 @@ elements = new DataSet(ElementsEnum()); nodes = new DataSet(NodesEnum()); + vertices = new DataSet(VerticesEnum()); materials = new DataSet(MaterialsEnum()); @@ -329,7 +331,12 @@ #endif + /*create vertex: */ + vertex_id=i+1; + vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i])); + + vertices->AddObject(vertex); + + node_id=i+1; //matlab indexing - - #ifdef _PARALLEL_ @@ -344,14 +351,10 @@ #endif - node_x[0]=iomodel->x[i]; - node_x[1]=iomodel->y[i]; - node_x[2]=iomodel->z[i]; - node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]); - - node_onbed=(int)iomodel->gridonbed[i]; - node_onsurface=(int)iomodel->gridonsurface[i]; - node_onshelf=(int)iomodel->gridoniceshelf[i]; - node_onsheet=(int)iomodel->gridonicesheet[i]; - + node_properties.SetProperties( + (int)iomodel->gridonbed[i], + (int)iomodel->gridonsurface[i], + (int)iomodel->gridoniceshelf[i], + (int)iomodel->gridonicesheet[i]); + if (isnan(iomodel->uppernodes[i])){ node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves. @@ -362,5 +365,5 @@ /*Create node using its constructor: */ - node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet); + node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties); /*Add node to nodes dataset: */ @@ -376,4 +379,5 @@ elements->Presort(); nodes->Presort(); + vertices->Presort(); materials->Presort(); @@ -409,4 +413,5 @@ *pelements=elements; *pnodes=nodes; + *pvertices=vertices; *pmaterials=materials; Index: /issm/trunk/src/c/ModelProcessorx/IoModel.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/IoModel.cpp (revision 3416) +++ /issm/trunk/src/c/ModelProcessorx/IoModel.cpp (revision 3417) @@ -46,5 +46,5 @@ iomodel->qmu_npart=0; iomodel->numberofelements=0; - iomodel->numberofnodes=0; + iomodel->numberofvertices=0; iomodel->x=NULL; iomodel->y=NULL; @@ -52,5 +52,5 @@ iomodel->elements=NULL; iomodel->elements_type=NULL; - iomodel->numberofnodes2d=0; + iomodel->numberofvertices2d=0; iomodel->elements2d=NULL; iomodel->deadgrids=NULL; @@ -180,9 +180,10 @@ iomodel->isstokes=0; - - iomodel->epart=NULL; - iomodel->npart=NULL; - iomodel->my_grids=NULL; - iomodel->my_bordergrids=NULL; + /*exterior data: */ + iomodel->my_elements=NULL; + iomodel->my_vertices=NULL; + iomodel->my_nodes=NULL; + iomodel->my_bordervertices=NULL; + iomodel->penaltypartitioning=NULL; return iomodel; @@ -276,9 +277,11 @@ xfree((void**)&iomodel->control_type); - xfree((void**)&iomodel->epart); - xfree((void**)&iomodel->npart); - xfree((void**)&iomodel->my_grids); - xfree((void**)&iomodel->my_bordergrids); - + /*exterior data: */ + xfree((void**)&iomodel->my_elements); + xfree((void**)&iomodel->my_vertices); + xfree((void**)&iomodel->my_nodes); + xfree((void**)&iomodel->my_bordervertices); + xfree((void**)&iomodel->penaltypartitioning); + /*!Delete entire structure: */ xfree((void**)piomodel); @@ -309,6 +312,6 @@ IoModelFetchData(&iomodel->control_analysis,iomodel_handle,"control_analysis"); IoModelFetchData(&iomodel->meshtype,iomodel_handle,"type"); - /*!Get numberofelements and numberofnodes: */ - IoModelFetchData(&iomodel->numberofnodes,iomodel_handle,"numberofgrids"); + /*!Get numberofelements and numberofvertices: */ + IoModelFetchData(&iomodel->numberofvertices,iomodel_handle,"numberofgrids"); IoModelFetchData(&iomodel->numberofelements,iomodel_handle,"numberofelements"); /*!In case we are running 3d, we are going to need the collapsed and non-collapsed 2d meshes, from which the 3d mesh was extruded: */ @@ -317,5 +320,5 @@ /*!Deal with 2d mesh: */ IoModelFetchData(&iomodel->numberofelements2d,iomodel_handle,"numberofelements2d"); - IoModelFetchData(&iomodel->numberofnodes2d,iomodel_handle,"numberofgrids2d"); + IoModelFetchData(&iomodel->numberofvertices2d,iomodel_handle,"numberofgrids2d"); IoModelFetchData(&iomodel->numlayers,iomodel_handle,"numlayers"); } Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 3416) +++ /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 3417) @@ -28,5 +28,5 @@ /*2d mesh: */ int numberofelements; - int numberofnodes; + int numberofvertices; double* x; double* y; @@ -36,5 +36,5 @@ /*3d mesh: */ - int numberofnodes2d; + int numberofvertices2d; int numberofelements2d; double* elements2d; @@ -176,9 +176,9 @@ int numoutput; - /*exterior data: */ - int* epart; /*!element partitioning.*/ - int* npart; /*!node partitioning.*/ - int* my_grids; /*! grids that belong to this cpu*/ - double* my_bordergrids; /*! grids that belong to this cpu, and some other cpu also*/ + /*exterior partitioning data, to be carried around: */ + bool* my_elements; + bool* my_vertices; + bool* my_nodes; + double* my_bordervertices; int* penaltypartitioning; @@ -197,5 +197,5 @@ /*Creation of fem datasets: general drivers*/ - void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateParameters(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); @@ -204,5 +204,5 @@ /*diagnostic horizontal*/ - void CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); @@ -210,20 +210,20 @@ /*diagnostic vertical*/ - void CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsDiagnosticVert(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsDiagnosticVert(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); /*diagnostic hutter*/ - void CreateElementsNodesAndMaterialsDiagnosticHutter(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateElementsNodesAndMaterialsDiagnosticHutter(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsDiagnosticHutter(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsDiagnosticHutter(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); /*diagnostic stokes*/ - void CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsDiagnosticStokes(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsDiagnosticStokes(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); /*slope compute*/ - void CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsSlopeCompute(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsSlopeCompute(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); @@ -233,5 +233,5 @@ /*thermal:*/ - void CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsThermal(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsThermal(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); @@ -239,5 +239,5 @@ /*melting:*/ - void CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsMelting(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsMelting(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); @@ -245,5 +245,5 @@ /*prognostic:*/ - void CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsPrognostic(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsPrognostic(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); @@ -251,5 +251,5 @@ /*prognostic2:*/ - void CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsPrognostic2(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsPrognostic2(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); @@ -257,5 +257,5 @@ /*balancedthickness:*/ - void CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsBalancedthickness(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); @@ -263,5 +263,5 @@ /*balancedvelocities:*/ - void CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsBalancedvelocities(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateLoadsBalancedvelocities(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); @@ -271,3 +271,7 @@ void CreateParametersQmu(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + + /*partitioning: */ + void Partitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, double** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle); + #endif /* _IOMODEL_H */ Index: /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp (revision 3416) +++ /issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp (revision 3417) @@ -13,5 +13,5 @@ -void CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ +void CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ @@ -26,8 +26,10 @@ DataSet* elements = NULL; DataSet* nodes = NULL; + DataSet* vertices = NULL; DataSet* materials = NULL; /*Objects: */ Node* node = NULL; + Vertex* vertex = NULL; Penta* penta = NULL; Matice* matice = NULL; @@ -88,7 +90,6 @@ /* node constructor input: */ int node_id; + int vertex_id; int node_partitionborder=0; - double node_x[3]; - double node_sigma; int node_onbed; int node_onsurface; @@ -113,4 +114,5 @@ elements = new DataSet(ElementsEnum()); nodes = new DataSet(NodesEnum()); + vertices = new DataSet(VerticesEnum()); materials = new DataSet(MaterialsEnum()); @@ -371,7 +373,11 @@ #endif + /*create vertex: */ + vertex_id=i+1; + vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i])); + vertices->AddObject(vertex); + + node_id=i+1; //matlab indexing - - #ifdef _PARALLEL_ @@ -386,13 +392,9 @@ #endif - node_x[0]=iomodel->x[i]; - node_x[1]=iomodel->y[i]; - node_x[2]=iomodel->z[i]; - node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]); - - node_onbed=(int)iomodel->gridonbed[i]; - node_onsurface=(int)iomodel->gridonsurface[i]; - node_onshelf=(int)iomodel->gridoniceshelf[i]; - node_onsheet=(int)iomodel->gridonicesheet[i]; + node_properties.SetProperties( + (int)iomodel->gridonbed[i], + (int)iomodel->gridonsurface[i], + (int)iomodel->gridoniceshelf[i], + (int)iomodel->gridonicesheet[i]); if (strcmp(iomodel->meshtype,"3d")==0){ @@ -410,5 +412,5 @@ /*Create node using its constructor: */ - node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet); + node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties); /*set single point constraints.: */ @@ -430,4 +432,5 @@ elements->Presort(); nodes->Presort(); + vertices->Presort(); materials->Presort(); @@ -462,4 +465,5 @@ *pelements=elements; *pnodes=nodes; + *pvertices=vertices; *pmaterials=materials; } Index: /issm/trunk/src/c/ModelProcessorx/Partitioning.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Partitioning.cpp (revision 3417) +++ /issm/trunk/src/c/ModelProcessorx/Partitioning.cpp (revision 3417) @@ -0,0 +1,172 @@ +/*!\file: Partitioning.cpp + * \brief: partition elements and nodes and vertices + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include +#include "./IoModel.h" + +void DiscontinuousGalerkinPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle); +void ContinuousGalerkinPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle); + +void Partitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){ + + if (iomodel->analysis_type==Prognostic2AnalysisEnum()) + DiscontinuousGalerkinPartitioning(pmy_elements, pmy_vertices, pmy_nodes, pmy_bordervertices, iomodel, iomodel_handle); + else + ContinuousGalerkinPartitioning(pmy_elements, pmy_vertices, pmy_nodes, pmy_bordervertices, iomodel, iomodel_handle); +} + +void ContinuousGalerkinPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){ + + /*as many nodes as there are vertices */ + + extern int my_rank; + extern int num_procs; + + /*output: */ + bool* my_elements=NULL; + bool* my_vertices=NULL; + bool* my_nodes=NULL; + bool* my_bordervertices=NULL; + + /*intermediary: */ + int* epart=NULL; //element partitioning. + int* npart=NULL; //node partitioning. + int elements_width; //size of elements + Vec bordervertices=NULL; + double* serial_bordervertices=NULL; + + /*Number of vertices per elements, needed to correctly retrieve data: */ + if(strcmp(iomodel->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(iomodel->meshtype,"2d")==0){ + /*load elements: */ + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + } + else{ + /*load elements2d: */ + IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d"); + } + + MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofvertices,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofvertices2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs); + + /*Free elements and elements2d: */ + xfree((void**)&iomodel->elements); + xfree((void**)&iomodel->elements2d); + + #else + /*In serial mode, epart is full of 0: all elements belong to cpu 0: */ + epart=(int*)xcalloc(md.numberofelements,sizeof(int)); + #endif + + /*Deal with rifts, they have to be included into one partition only, not several: */ + IoModelFetchData(&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo"); + + for(i=0;inumrifts;i++){ + el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2)-1; //matlab indexing to c indexing + el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+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**)&iomodel->riftinfo); + + /*Used later on: */ + my_vertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool)); + my_elements=(boll*)xcalloc(iomodel->numberofelements,sizeof(bool)); + + /*Start figuring out, out of the partition, which elements belong to this cpu: */ + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + for (i=0;inumberofelements;i++){ + + /*!All elements have been partitioned above, only deal with elements for this cpu: */ + if(my_rank==epart[i]){ + + my_elements[i]=1; + + /*Now that we are here, we can also start building the list of vertices belonging to this cpu 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 vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices + will hold which vertices belong to this partition*/ + my_vertices[(int)*(iomodel->elements+elements_width*i+0)-1]=1; + my_vertices[(int)*(iomodel->elements+elements_width*i+1)-1]=1; + my_vertices[(int)*(iomodel->elements+elements_width*i+2)-1]=1; + + if(elements_width==6){ + my_vertices[(int)*(iomodel->elements+elements_width*i+3)-1]=1; + my_vertices[(int)*(iomodel->elements+elements_width*i+4)-1]=1; + my_vertices[(int)*(iomodel->elements+elements_width*i+5)-1]=1; + } + } + }//for (i=0;ielements); + + #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:*/ + bordervertices=NewVec(iomodel->numberofvertices); + + for (i=0;inumberofvertices;i++){ + if(my_vertices[i])VecSetValue(bordervertices,i,1,ADD_VALUES); + } + VecAssemblyBegin(bordervertices); + VecAssemblyEnd(bordervertices); + + VecToMPISerial(&serial_bordervertices,bordervertices); + + /*now go through serial_bordervertices, and booleanize it: */ + my_bordervertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool)); + for(i=0;inumberofvertices;i++){ + if(serial_bordervertices>1)my_bordervertices[i]=1; + } + #else + /*No border vertices: */ + my_bordervertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool)); + #endif + + /*Deal with my_nodes: */ + my_nodes=(bool*)xmalloc(iomodel->numberofnodes*sizeof(bool)); + memcpy(my_nodes,my_vertices,iomodel->numberofnodes*sizeof(bool)); + + /*Free ressources:*/ + xfree((void**)&npart); + xfree((void**)&epart); + VecFree(&bordervertices); + + /*Assign output pointers:*/ + *pmy_elements=my_elements; + *pmy_vertices=my_vertices; + *pmy_nodes=my_nodes; + *pmy_bordervertices=my_bordervertices; +} + + +void DiscontinuousGalerkinPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){ + + /*each element has it own nodes (as many as vertices) + additional nodes from neighbouring elements for each edge. This yields to a very different partition for + * the nodes and the vertices. The vertices are similar to continuous galerkin, but the nodes partitioning involves edges, which mess up sorting of + * ids. */ + + /*output: */ + bool* my_elements=NULL; + bool* my_vertices=NULL; + bool* my_nodes=NULL; + bool* my_bordervertices=NULL; + + /*Mathieu's code: */ + + /*Assign output pointers:*/ + *pmy_elements=my_elements; + *pmy_vertices=my_vertices; + *pmy_nodes=my_nodes; + *pmy_bordervertices=my_bordervertices; +} Index: /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp (revision 3416) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp (revision 3417) @@ -14,5 +14,5 @@ -void CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ +void CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ @@ -27,8 +27,10 @@ DataSet* elements = NULL; DataSet* nodes = NULL; + DataSet* vertices = NULL; DataSet* materials = NULL; /*Objects: */ Node* node = NULL; + Vertex* vertex = NULL; Tria* tria = NULL; Penta* penta = NULL; @@ -97,7 +99,6 @@ /* node constructor input: */ int node_id; + int vertex_id; int node_partitionborder=0; - double node_x[3]; - double node_sigma; int node_onbed; int node_onsurface; @@ -122,4 +123,5 @@ elements = new DataSet(ElementsEnum()); nodes = new DataSet(NodesEnum()); + vertices = new DataSet(VerticesEnum()); materials = new DataSet(MaterialsEnum()); @@ -436,7 +438,11 @@ #endif + /*create vertex: */ + vertex_id=i+1; + vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i])); + vertices->AddObject(vertex); + + node_id=i+1; //matlab indexing - - #ifdef _PARALLEL_ @@ -451,14 +457,11 @@ #endif - node_x[0]=iomodel->x[i]; - node_x[1]=iomodel->y[i]; - node_x[2]=iomodel->z[i]; - node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]); - - node_onbed=(int)iomodel->gridonbed[i]; - node_onsurface=(int)iomodel->gridonsurface[i]; - node_onshelf=(int)iomodel->gridoniceshelf[i]; - node_onsheet=(int)iomodel->gridonicesheet[i]; - + node_properties.SetProperties( + (int)iomodel->gridonbed[i], + (int)iomodel->gridonsurface[i], + (int)iomodel->gridoniceshelf[i], + (int)iomodel->gridonicesheet[i]); + + if (strcmp(iomodel->meshtype,"3d")==0){ if (isnan(iomodel->uppernodes[i])){ @@ -475,5 +478,5 @@ /*Create node using its constructor: */ - node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet); + node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties); /*set single point constraints.: */ @@ -498,4 +501,5 @@ elements->Presort(); nodes->Presort(); + vertices->Presort(); materials->Presort(); @@ -531,4 +535,5 @@ *pelements=elements; *pnodes=nodes; + *pvertices=vertices; *pmaterials=materials; Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3416) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3417) @@ -13,5 +13,5 @@ #include "../IoModel.h" -void CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ +void CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ int i,j,k,n; @@ -22,8 +22,10 @@ DataSet* elements = NULL; DataSet* nodes = NULL; + DataSet* vertices = NULL; DataSet* materials = NULL; /*Objects: */ Node* node = NULL; + Vertex* vertex = NULL; Tria* tria = NULL; Penta* penta = NULL; @@ -96,7 +98,7 @@ /* node constructor input: */ int node_id; + int vertex_id; + int vertex_partitionborder=0; int node_partitionborder=0; - double node_x[3]; - double node_sigma; int node_onbed; int node_onsurface; @@ -125,4 +127,5 @@ elements = new DataSet(ElementsEnum()); nodes = new DataSet(NodesEnum()); + vertices = new DataSet(VerticesEnum()); materials = new DataSet(MaterialsEnum()); @@ -349,4 +352,22 @@ DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); + /*Create vertices: */ + for (i=0;inumberofnodes;i++){ + #ifdef _PARALLEL_ + /*keep only this partition's nodes:*/ + if((my_grids[i]==1)){ + #endif + + /*create vertex: */ + vertex_id=i+1; + vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i])); + + vertices->AddObject(vertex); + + #ifdef _PARALLEL_ + } //if((my_grids[i]==1)) + #endif + } + /*Build Nodes dataset -> 3 for each element*/ for (i=0;inumberofelements;i++){ @@ -374,12 +395,11 @@ node_partitionborder=0; #endif - node_x[0]=iomodel->x[k]; - node_x[1]=iomodel->y[k]; - node_x[2]=iomodel->z[k]; - node_sigma=(iomodel->z[k]-iomodel->bed[k])/(iomodel->thickness[k]); - node_onbed=(int)iomodel->gridonbed[k]; - node_onsurface=(int)iomodel->gridonsurface[k]; - node_onshelf=(int)iomodel->gridoniceshelf[k]; - node_onsheet=(int)iomodel->gridonicesheet[k]; + + node_properties.SetProperties( + (int)iomodel->gridonbed[k], + (int)iomodel->gridonsurface[k], + (int)iomodel->gridoniceshelf[k], + (int)iomodel->gridonicesheet[k]); + if (strcmp(iomodel->meshtype,"3d")==0){ if (isnan(iomodel->uppernodes[k])){ @@ -396,5 +416,5 @@ /*Create node using its constructor: */ - node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet); + node=new Node(node_id, k, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties); /*Add node to nodes dataset: */ @@ -410,4 +430,5 @@ elements->Presort(); nodes->Presort(); + vertices->Presort(); materials->Presort(); @@ -443,4 +464,5 @@ *pelements=elements; *pnodes=nodes; + *pvertices=vertices; *pmaterials=materials; } Index: /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 3416) +++ /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 3417) @@ -13,5 +13,5 @@ #include "../IoModel.h" -void CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ +void CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ @@ -26,8 +26,10 @@ DataSet* elements = NULL; DataSet* nodes = NULL; + DataSet* vertices = NULL; DataSet* materials = NULL; /*Objects: */ Node* node = NULL; + Vertex* vertex = NULL; Tria* tria = NULL; Penta* penta = NULL; @@ -68,7 +70,6 @@ /* node constructor input: */ int node_id; + int vertex_id; int node_partitionborder=0; - double node_x[3]; - double node_sigma; int node_onbed; int node_onsurface; @@ -93,4 +94,5 @@ elements = new DataSet(ElementsEnum()); nodes = new DataSet(NodesEnum()); + vertices = new DataSet(VerticesEnum()); materials = new DataSet(MaterialsEnum()); @@ -372,7 +374,11 @@ #endif + /*create vertex: */ + vertex_id=i+1; + vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i])); + vertices->AddObject(vertex); + + node_id=i+1; //matlab indexing - - #ifdef _PARALLEL_ @@ -387,13 +393,9 @@ #endif - node_x[0]=iomodel->x[i]; - node_x[1]=iomodel->y[i]; - node_x[2]=iomodel->z[i]; - node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]); - - node_onbed=(int)iomodel->gridonbed[i]; - node_onsurface=(int)iomodel->gridonsurface[i]; - node_onshelf=(int)iomodel->gridoniceshelf[i]; - node_onsheet=(int)iomodel->gridonicesheet[i]; + node_properties.SetProperties( + (int)iomodel->gridonbed[i], + (int)iomodel->gridonsurface[i], + (int)iomodel->gridoniceshelf[i], + (int)iomodel->gridonicesheet[i]); if (strcmp(iomodel->meshtype,"3d")==0){ @@ -411,5 +413,5 @@ /*Create node using its constructor: */ - node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet); + node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties); /*set single point constraints.: */ @@ -434,4 +436,5 @@ elements->Presort(); nodes->Presort(); + vertices->Presort(); materials->Presort(); @@ -467,4 +470,5 @@ *pelements=elements; *pnodes=nodes; + *pvertices=vertices; *pmaterials=materials; Index: /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp (revision 3416) +++ /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp (revision 3417) @@ -12,5 +12,5 @@ #include "../IoModel.h" -void CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ +void CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ @@ -25,8 +25,10 @@ DataSet* elements = NULL; DataSet* nodes = NULL; + DataSet* vertices = NULL; DataSet* materials = NULL; /*Objects: */ Node* node = NULL; + Vertex* vertex = NULL; Penta* penta = NULL; Matice* matice = NULL; @@ -86,7 +88,6 @@ /* node constructor input: */ int node_id; + int vertex_id; int node_partitionborder=0; - double node_x[3]; - double node_sigma; int node_onbed; int node_onsurface; @@ -111,4 +112,5 @@ elements = new DataSet(ElementsEnum()); nodes = new DataSet(NodesEnum()); + vertices = new DataSet(VerticesEnum()); materials = new DataSet(MaterialsEnum()); @@ -369,4 +371,10 @@ #endif + /*create vertex: */ + vertex_id=i+1; + vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i])); + vertices->AddObject(vertex); + + node_id=i+1; //matlab indexing @@ -382,13 +390,10 @@ #endif - node_x[0]=iomodel->x[i]; - node_x[1]=iomodel->y[i]; - node_x[2]=iomodel->z[i]; - node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]); - - node_onbed=(int)iomodel->gridonbed[i]; - node_onsurface=(int)iomodel->gridonsurface[i]; - node_onshelf=(int)iomodel->gridoniceshelf[i]; - node_onsheet=(int)iomodel->gridonicesheet[i]; + node_properties.SetProperties( + (int)iomodel->gridonbed[i], + (int)iomodel->gridonsurface[i], + (int)iomodel->gridoniceshelf[i], + (int)iomodel->gridonicesheet[i]); + if (strcmp(iomodel->meshtype,"3d")==0){ @@ -406,5 +411,5 @@ /*Create node using its constructor: */ - node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet); + node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties); /*Add node to nodes dataset: */ @@ -420,4 +425,5 @@ elements->Presort(); nodes->Presort(); + vertices->Presort(); materials->Presort(); @@ -440,4 +446,5 @@ *pelements=elements; *pnodes=nodes; + *pvertices=vertices; *pmaterials=materials; Index: /issm/trunk/src/c/objects/DofIndexing.cpp =================================================================== --- /issm/trunk/src/c/objects/DofIndexing.cpp (revision 3417) +++ /issm/trunk/src/c/objects/DofIndexing.cpp (revision 3417) @@ -0,0 +1,183 @@ +/*!\file DofIndexing.c + * \brief: implementation of the DofIndexing object + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "stdio.h" +#include "./DofIndexing.h" +#include +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../include/typedefs.h" +#include "../include/macros.h" + +/*Object constructors and destructor*/ +/*FUNCTION DofIndexing default constructor {{{1*/ +DofIndexing::DofIndexing(){ + + int i; + this->numberofdofs=UNDEF; + this->partitionborder=UNDEF; + this->clone=UNDEF; + + for (i=0;im_set[i]=0; + this->n_set[i]=1; + this->f_set[i]=1; + this->s_set[i]=0; + this->doflist[i]=UNDEF; + } +} +/*}}}*/ +/*FUNCTION DofIndexing constructor {{{1*/ +DofIndexing::DofIndexing(int in_numberofdofs, int in_partitionborder){ + + int i; + this->numberofdofs=in_numberofdofs; + this->partitionborder=in_partitionborder; + this->clone=UNDEF; + + for (i=0;im_set[i]=0; + this->n_set[i]=1; + this->f_set[i]=1; + this->s_set[i]=0; + this->doflist[i]=UNDEF; + } +} +/*}}}*/ +/*FUNCTION DofIndexing copy constructor{{{1*/ +DofIndexing::DofIndexing(DofIndexing* in){ //copy constructor + + int i; + this->numberofdofs=in->numberofdofs; + this->partitionborder=in->partitionborder; + this->clone=in->clone; + + for(i=0;im_set[i]=in->m_set[i]; + this->n_set[i]=in->n_set[i]; + this->f_set[i]=in->f_set[i]; + this->s_set[i]=in->s_set[i]; + this->doflist[i]=in->doflist[i]; + } +} +/*}}}*/ +/*FUNCTION DofIndexing destructor {{{1*/ +DofIndexing::~DofIndexing(){ //destructor +} +/*}}}*/ + +/*Object management: */ +/*FUNCTION DofIndexing Marshall{{{1*/ +void DofIndexing::Marshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int enum_type=0; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*get enum type of DofIndexing: */ + enum_type=DofIndexingEnum(); + + /*marshall enum: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + + /*marshall DofIndexing data: */ + memcpy(marshalled_dataset,&numberofdofs,sizeof(numberofdofs));marshalled_dataset+=sizeof(numberofdofs); + memcpy(marshalled_dataset,&partitionborder,sizeof(partitionborder));marshalled_dataset+=sizeof(partitionborder); + memcpy(marshalled_dataset,&clone,sizeof(clone));marshalled_dataset+=sizeof(clone); + memcpy(marshalled_dataset,&m_set,sizeof(m_set));marshalled_dataset+=sizeof(m_set); + memcpy(marshalled_dataset,&n_set,sizeof(n_set));marshalled_dataset+=sizeof(n_set); + memcpy(marshalled_dataset,&f_set,sizeof(f_set));marshalled_dataset+=sizeof(f_set); + memcpy(marshalled_dataset,&s_set,sizeof(s_set));marshalled_dataset+=sizeof(s_set); + memcpy(marshalled_dataset,&doflist,sizeof(doflist));marshalled_dataset+=sizeof(doflist); + + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ +/*FUNCTION DofIndexing MarshallSize{{{1*/ +int DofIndexing::MarshallSize(){ + + return sizeof(numberofdofs)+ + sizeof(partitionborder)+ + sizeof(clone)+ + sizeof(m_set)+ + sizeof(n_set)+ + sizeof(f_set)+ + sizeof(s_set)+ + sizeof(doflist)+ + sizeof(int); //sizeof(int) for enum type +} +/*}}}*/ +/*FUNCTION DofIndexing Demarshall{{{1*/ +void DofIndexing::Demarshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*this time, no need to get enum type, the pointer directly points to the beginning of the + *object data (thanks to DataSet::Demarshall):*/ + + memcpy(&numberofdofs,marshalled_dataset,sizeof(numberofdofs));marshalled_dataset+=sizeof(numberofdofs); + memcpy(&partitionborder,marshalled_dataset,sizeof(partitionborder));marshalled_dataset+=sizeof(partitionborder); + memcpy(&clone,marshalled_dataset,sizeof(clone));marshalled_dataset+=sizeof(clone); + memcpy(&m_set,marshalled_dataset,sizeof(m_set));marshalled_dataset+=sizeof(m_set); + memcpy(&n_set,marshalled_dataset,sizeof(n_set));marshalled_dataset+=sizeof(n_set); + memcpy(&f_set,marshalled_dataset,sizeof(f_set));marshalled_dataset+=sizeof(f_set); + memcpy(&s_set,marshalled_dataset,sizeof(s_set));marshalled_dataset+=sizeof(s_set); + memcpy(&doflist,marshalled_dataset,sizeof(doflist));marshalled_dataset+=sizeof(doflist); + + /*return: */ + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ +/*FUNCTION DofIndexing Echo{{{1*/ +void DofIndexing::Echo(void){ + + int i; + + printf("DofIndexing:\n"); + printf(" numberofdofs: %i\n",numberofdofs); + printf(" partitionborder: %i\n",partitionborder); + printf(" clone: %i\n",clone); +} +/*}}}*/ +/*FUNCTION DofIndexing DeepEcho{{{1*/ +void DofIndexing::DeepEcho(void){ + + int i; + + printf("DofIndexing:\n"); + printf(" numberofdofs: %i\n",numberofdofs); + printf(" partitionborder: %i\n",partitionborder); + printf(" clone: %i\n",clone); + + printf(" set membership: m,n,f,s sets \n"); + for(i=0;iMAXDOFSPERNODE)break; + printf(" dof %i: %i %i %i %i\n",i,m_set[i],n_set[i],f_set[i],s_set[i]); + } + printf("\n"); + + printf(" doflist:|"); + for(i=0;iMAXDOFSPERNODE)break; + printf("%i|",doflist[i]); + } + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/DofIndexing.h =================================================================== --- /issm/trunk/src/c/objects/DofIndexing.h (revision 3417) +++ /issm/trunk/src/c/objects/DofIndexing.h (revision 3417) @@ -0,0 +1,45 @@ +/*!\file: DofIndexing.h + * \brief prototype for DofIndexing.h + */ + +#ifndef _DOFINDEXING_H_ +#define _DOFINDEXING_H_ + +#define MAXDOFSPERNODE 4 + +class DofIndexing{ + + public: + + int numberofdofs; //number of dofs for a node. + + /*partitioning: */ + int partitionborder; /*! during parallel partitioning, does this grid belong to a partition border?*/ + int clone; /*!this nodes is one the partition border, and is cloned*/ + + /*boundary conditions sets: */ + int m_set[MAXDOFSPERNODE]; + int n_set[MAXDOFSPERNODE]; + int f_set[MAXDOFSPERNODE]; + int s_set[MAXDOFSPERNODE]; + + /*list of degrees of freedom: */ + int doflist[MAXDOFSPERNODE]; //dof list on which we solve + + DofIndexing(); + DofIndexing(int numberofdofs, int partitionborder); + DofIndexing(DofIndexing* properties); + ~DofIndexing(); + + void Marshall(char** pmarshalled_dataset); + int MarshallSize(); + void Demarshall(char** pmarshalled_dataset); + void copy(DofIndexing* properties); + void Echo(void); + void DeepEcho(void); + DofIndexing* Spawn(int* indices, int numindices); + +}; +#endif //ifndef _DOFINDEXING_H_ + + Index: /issm/trunk/src/c/objects/ElementProperties.cpp =================================================================== --- /issm/trunk/src/c/objects/ElementProperties.cpp (revision 3416) +++ /issm/trunk/src/c/objects/ElementProperties.cpp (revision 3417) @@ -96,7 +96,26 @@ } /*}}}*/ -/*FUNCTION ElementProperties copy constructor{{{1*/ -ElementProperties::ElementProperties(ElementProperties* prop){ //copy constructor - +/*FUNCTION ElementProperties constructor {{{1*/ + +ElementProperties::ElementProperties(int elementproperties_numnodes, double* elementproperties_h,double* elementproperties_s,double* elementproperties_b, + double* elementproperties_k,double* elementproperties_melting,double* elementproperties_accumulation, + double* elementproperties_geothermalflux, int elementproperties_friction_type,double elementproperties_p, + double elementproperties_q, int elementproperties_shelf, int elementproperties_onbed, bool elementproperties_onwater, + int elementproperties_onsurface, int elementproperties_collapse, int elementproperties_thermal_steadystate){ + + this->Init(elementproperties_numnodes, elementproperties_h,elementproperties_s,elementproperties_b, + elementproperties_k,elementproperties_melting,elementproperties_accumulation, + elementproperties_geothermalflux, elementproperties_friction_type,elementproperties_p, + elementproperties_q, elementproperties_shelf, elementproperties_onbed, elementproperties_onwater, + elementproperties_onsurface, elementproperties_collapse, elementproperties_thermal_steadystate); +} +/*}}}*/ +/*FUNCTION ElementProperties Initialize propreties, used by constructor{{{1*/ +ElementProperties::Init(int elementproperties_numnodes, double* elementproperties_h,double* elementproperties_s,double* elementproperties_b, + double* elementproperties_k,double* elementproperties_melting,double* elementproperties_accumulation, + double* elementproperties_geothermalflux, int elementproperties_friction_type,double elementproperties_p, + double elementproperties_q, int elementproperties_shelf, int elementproperties_onbed, bool elementproperties_onwater, + int elementproperties_onsurface, int elementproperties_collapse, int elementproperties_thermal_steadystate){ + int i; Index: /issm/trunk/src/c/objects/ElementProperties.h =================================================================== --- /issm/trunk/src/c/objects/ElementProperties.h (revision 3416) +++ /issm/trunk/src/c/objects/ElementProperties.h (revision 3417) @@ -33,4 +33,5 @@ ElementProperties(); ElementProperties(int numnodes, double* h, double* s, double* b, double* k, double* melting, double* accumulation, double* geothermalflux, int friction_type, double p, double q, int shelf, int onbed, bool onwater, int onsurface, int collapse, int thermal_steadystate); + Init(int numnodes, double* h, double* s, double* b, double* k, double* melting, double* accumulation, double* geothermalflux, int friction_type, double p, double q, int shelf, int onbed, bool onwater, int onsurface, int collapse, int thermal_steadystate); ElementProperties(ElementProperties* properties); ~ElementProperties(); Index: /issm/trunk/src/c/objects/FemModel.cpp =================================================================== --- /issm/trunk/src/c/objects/FemModel.cpp (revision 3416) +++ /issm/trunk/src/c/objects/FemModel.cpp (revision 3417) @@ -9,13 +9,15 @@ #endif -#include "./FemModel.h" #include "stdio.h" #include "../shared/shared.h" #include "../include/macros.h" - +#include "./FemModel.h" + +/*Object constructors, destructors: {{{1:*/ FemModel::FemModel(){ elements=NULL; nodes=NULL; + vertices=NULL; constraints=NULL; loads=NULL; @@ -34,5 +36,5 @@ } -FemModel::FemModel(DataSet* femmodel_elements,DataSet* femmodel_nodes,DataSet* femmodel_constraints,DataSet* femmodel_loads, +FemModel::FemModel(DataSet* femmodel_elements,DataSet* femmodel_nodes,DataSet* femmodel_vertices, DataSet* femmodel_constraints,DataSet* femmodel_loads, DataSet* femmodel_materials,DataSet* femmodel_parameters, DofVec* femmodel_partition,DofVec* femmodel_tpartition,DofVec* femmodel_yg, Mat femmodel_Rmg,Mat femmodel_Gmn,NodeSets* femmodel_nodesets,Vec femmodel_ys,Vec femmodel_ys0){ @@ -41,4 +43,5 @@ elements=femmodel_elements; nodes=femmodel_nodes; + vertices=femmodel_vertices; constraints=femmodel_constraints; loads=femmodel_loads; @@ -61,4 +64,5 @@ delete elements; delete nodes; + delete vertices; delete loads; delete constraints; @@ -76,6 +80,6 @@ } - - +/*}}}*/ +/* Object management: {{{1*/ void FemModel::Echo(void){ @@ -83,4 +87,5 @@ printf(" elements: %p\n",elements); printf(" nodes: %p\n",nodes); + printf(" vertices: %p\n",vertices); printf(" loads: %p\n",loads); printf(" materials: %p\n",materials); @@ -105,4 +110,6 @@ printf(" nodes: \n"); nodes->Echo(); + printf(" vertices: \n"); + vertices->Echo(); printf(" loads: \n"); nodes->Echo(); @@ -204,4 +211,5 @@ DataSet* FemModel::get_elements(void){return elements;} DataSet* FemModel::get_nodes(void){return nodes ;} +DataSet* FemModel::get_vertices(void){return vertices ;} DataSet* FemModel::get_constraints(void){return constraints ;} DataSet* FemModel::get_loads(void){return loads;} @@ -216,2 +224,3 @@ Vec FemModel::get_ys0(void){return ys0;} Mat FemModel::get_Gmn(void){return Gmn;} +/*}}}*/ Index: /issm/trunk/src/c/objects/FemModel.h =================================================================== --- /issm/trunk/src/c/objects/FemModel.h (revision 3416) +++ /issm/trunk/src/c/objects/FemModel.h (revision 3417) @@ -6,11 +6,9 @@ #define FEMMODEL_H_ +#include "./Object.h" +#include "../DataSet/DataSet.h" +#include "./DofVec.h" #include "../toolkits/toolkits.h" -#include "../DataSet/DataSet.h" #include "../parallel/parallel.h" -#include "./DofVec.h" - -class DataSet; -struct OptArgs; class FemModel: public Object{ @@ -22,4 +20,5 @@ DataSet* elements; DataSet* nodes; + DataSet* vertices; DataSet* constraints; DataSet* loads; @@ -27,7 +26,8 @@ DataSet* parameters; - DofVec* partition; - DofVec* tpartition; - DofVec* yg; + DofVec* partition; + DofVec* tpartition; + DofVec* yg; + Mat Rmg; NodeSets* nodesets; @@ -38,5 +38,5 @@ FemModel(); ~FemModel(); - FemModel(DataSet* elements,DataSet* nodes,DataSet* constraints,DataSet* loads,DataSet* materials,DataSet* parameters, + FemModel(DataSet* elements,DataSet* nodes,DataSet* vertices, DataSet* constraints,DataSet* loads,DataSet* materials,DataSet* parameters, DofVec* partition,DofVec* tpartition,DofVec* yg,Mat Rmg,Mat Gmn,NodeSets* nodesets,Vec ys,Vec ys0); @@ -62,4 +62,5 @@ DataSet* get_elements(void); DataSet* get_nodes(void); + DataSet* get_vertices(void){return vertices ;} DataSet* get_constraints(void); DataSet* get_loads(void); @@ -75,7 +76,4 @@ Mat get_Gmn(void); - - - }; Index: /issm/trunk/src/c/objects/Hook.cpp =================================================================== --- /issm/trunk/src/c/objects/Hook.cpp (revision 3416) +++ /issm/trunk/src/c/objects/Hook.cpp (revision 3417) @@ -12,4 +12,7 @@ #include "stdio.h" #include +#include "./Object.h" +#include "../DataSet/DataSet.h" +#include "./Hook.h" #include "../EnumDefinitions/EnumDefinitions.h" #include "./ParameterInputs.h" @@ -17,5 +20,4 @@ #include "../include/typedefs.h" #include "../include/macros.h" -#include "./Hook.h" @@ -31,8 +33,13 @@ /*}}}*/ /*FUNCTION Hook::Hook(int* ids, int num){{{1*/ -Hook::Hook(int* ids, int num){ - - int i; - this->num=num; +Hook::Hook(int* in_ids, int in_num){ + this->Init(in_ids,in_num); +} +/*}}}*/ +/*FUNCTION Hook::Init(int* ids, int num){{{1*/ +Hook::Init(int* in_ids, int in_num){ + + int i; + this->num=in_num; /*Allocate: */ @@ -43,5 +50,5 @@ /*Copy ids: */ for (i=0;inum;i++){ - this->ids[i]=ids[i]; + this->ids[i]=in_ids[i]; this->objects[i]=NULL; this->offsets[i]=0; Index: /issm/trunk/src/c/objects/Hook.h =================================================================== --- /issm/trunk/src/c/objects/Hook.h (revision 3416) +++ /issm/trunk/src/c/objects/Hook.h (revision 3417) @@ -10,10 +10,8 @@ #include "./Object.h" -#include "./../DataSet/DataSet.h" +#include "../DataSet/DataSet.h" #include "../toolkits/toolkits.h" class DataSet; -class Object; - class Hook{ @@ -29,4 +27,5 @@ Hook(); Hook(int* ids, int num); + Init(int* ids, int num); Hook(Object** objects, int* ids, int* offsets,int num); Hook(Hook* hook); Index: /issm/trunk/src/c/objects/Matice.cpp =================================================================== --- /issm/trunk/src/c/objects/Matice.cpp (revision 3416) +++ /issm/trunk/src/c/objects/Matice.cpp (revision 3417) @@ -17,15 +17,48 @@ /*Object constructors and destructor*/ +/*FUNCTION Matice::default constructor {{{1*/ +Matice::Matice(){ + return; +} +/*}}}*/ /*FUNCTION Matice::constructor {{{1*/ -Matice::Matice(){ - return; -} -/*}}}*/ -/*FUNCTION Matice::creation {{{1*/ -Matice::Matice(int matice_mid,double matice_B,double matice_n){ - mid=matice_mid; - B=matice_B; - n=matice_n; - return; +Matice::Matice(int in_mid,double in_B,double in_n){ + this->Init(in_mid,in_B,in_n); +} +/*}}}*/ +/*FUNCTION Matice::init {{{1*/ +Matice::Init(int in_mid,double in_B,double in_n){ + this->mid=in_mid; + this->B=in_B; + this->n=in_n; +} +/*}}}*/ +/*FUNCTION Matice::constructor from iomodel {{{1*/ +Matice::Matice(int i, IoModel* iomodel, int num_vertices){ + + int j; + + /*needed for Init routine:*/ + int matice_mid; + double matice_B; + double matice_n; + + /*intermediary: */ + double B_avg; + + /*id: */ + matice_mid=i+1; //same as element it is created for + + /*B and n: */ + B_avg=0; + for(j=0;jB+((int)*(iomodel->elements+num_vertices*i+j)-1)); + } + B_avg=B_avg/num_vertices; + + matice_B=B_avg; + matice_n=(double)*(iomodel->n+i); + + this->Init(matice_mid,matice_B,matice_n); } /*}}}*/ Index: /issm/trunk/src/c/objects/Matice.h =================================================================== --- /issm/trunk/src/c/objects/Matice.h (revision 3416) +++ /issm/trunk/src/c/objects/Matice.h (revision 3417) @@ -19,4 +19,6 @@ Matice(); Matice(int mid,double B,double n); + Matice(int i, IoModel* iomodel, int num_vertices); + Init(int mid,double B,double n); ~Matice(); Index: /issm/trunk/src/c/objects/Matpar.cpp =================================================================== --- /issm/trunk/src/c/objects/Matpar.cpp (revision 3416) +++ /issm/trunk/src/c/objects/Matpar.cpp (revision 3417) @@ -18,24 +18,61 @@ /*Object constructors and destructor*/ -/*FUNCTION Matpar::constructor {{{1*/ +/*FUNCTION Matpar::default constructor {{{1*/ Matpar::Matpar(){ return; } /*}}}1*/ -/*FUNCTION Matpar::creation {{{1*/ +/*FUNCTION Matpar::constructorr {{{1*/ Matpar::Matpar(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){ - - mid=matpar_mid; - rho_ice=matpar_rho_ice; - rho_water=matpar_rho_water; - heatcapacity=matpar_heatcapacity; - thermalconductivity=matpar_thermalconductivity; - latentheat=matpar_latentheat; - beta=matpar_beta; - meltingpoint=matpar_meltingpoint; - mixed_layer_capacity=matpar_mixed_layer_capacity; - thermal_exchange_velocity=matpar_thermal_exchange_velocity; - g=matpar_g; + this->Init(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); + +} +/*}}}1*/ +/*FUNCTION Matpar::constructor from iomodel{{{1*/ +Matpar::Matpar(IoModel* iomodel){ + + 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; + + matpar_mid=iomodel->numberofelements+1; //put it at the end of the materials + matpar_g=iomodel->g; + matpar_rho_ice=iomodel->rho_ice; + matpar_rho_water=iomodel->rho_water; + matpar_thermalconductivity=iomodel->thermalconductivity; + matpar_heatcapacity=iomodel->heatcapacity; + matpar_latentheat=iomodel->latentheat; + matpar_beta=iomodel->beta; + matpar_meltingpoint=iomodel->meltingpoint; + matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; + matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; + + this->Init(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); + +} +/*}}}1*/ +/*FUNCTION Matpar::Init {{{1*/ +Matpar::Init(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){ + + this->mid=matpar_mid; + this->rho_ice=matpar_rho_ice; + this->rho_water=matpar_rho_water; + this->heatcapacity=matpar_heatcapacity; + this->thermalconductivity=matpar_thermalconductivity; + this->latentheat=matpar_latentheat; + this->beta=matpar_beta; + this->meltingpoint=matpar_meltingpoint; + this->mixed_layer_capacity=matpar_mixed_layer_capacity; + this->thermal_exchange_velocity=matpar_thermal_exchange_velocity; + this->g=matpar_g; return; Index: /issm/trunk/src/c/objects/Matpar.h =================================================================== --- /issm/trunk/src/c/objects/Matpar.h (revision 3416) +++ /issm/trunk/src/c/objects/Matpar.h (revision 3417) @@ -28,4 +28,6 @@ Matpar(int mid, double rho_ice, double rho_water, double heatcapacity, double thermalconductivity, double latentheat, double beta, double meltingpoint, double mixed_layer_capacity, double thermal_exchange_velocity, double g); + Init(int mid, double rho_ice, double rho_water, double heatcapacity, double thermalconductivity, double latentheat, double beta, double meltingpoint, double mixed_layer_capacity, double thermal_exchange_velocity, double g); + Matpar(IoModel* iomodel); ~Matpar(); Index: /issm/trunk/src/c/objects/Model.cpp =================================================================== --- /issm/trunk/src/c/objects/Model.cpp (revision 3416) +++ /issm/trunk/src/c/objects/Model.cpp (revision 3417) @@ -43,4 +43,5 @@ DataSet* elements=NULL; DataSet* nodes=NULL; + DataSet* vertices=NULL; DataSet* constraints=NULL; DataSet* loads=NULL; @@ -71,8 +72,8 @@ _printf_(" create datasets:\n"); - CreateDataSets(&elements,&nodes,&materials,&constraints,&loads,¶meters,iomodel,MODEL); + CreateDataSets(&elements,&nodes,&vertices,&materials,&constraints,&loads,¶meters,iomodel,MODEL); _printf_(" create degrees of freedom: \n"); - Dofx( &partition,&tpartition,elements,nodes, parameters); + Dofx( &partition,&tpartition,elements,nodes, vertices,parameters); _printf_(" create single point constraints: \n"); @@ -92,5 +93,5 @@ _printf_(" configuring element and loads:\n"); - ConfigureObjectsx(elements, loads, nodes, materials,parameters); + ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters); _printf_(" process parameters:\n"); @@ -101,5 +102,5 @@ /*Use all the data created to create an femmodel: */ - femmodel=new FemModel(elements,nodes,constraints,loads,materials,parameters, + femmodel=new FemModel(elements, nodes, vertices, constraints,loads,materials,parameters, partition,tpartition,yg,Rmg,Gmn,nodesets,ys,ys0); @@ -116,4 +117,5 @@ DataSet* elements=NULL; DataSet* nodes=NULL; + DataSet* vertices=NULL; DataSet* constraints=NULL; DataSet* loads=NULL; @@ -148,8 +150,8 @@ _printf_(" create datasets:\n"); - CreateDataSets(&elements,&nodes,&materials,&constraints,&loads,¶meters,iomodel,MODEL); + CreateDataSets(&elements,&nodes,&vertices,&materials,&constraints,&loads,¶meters,iomodel,MODEL); _printf_(" create degrees of freedom: \n"); - Dofx( &partition,&tpartition,elements,nodes, parameters); + Dofx( &partition,&tpartition,elements,nodes, vertices,parameters); _printf_(" create single point constraints: \n"); @@ -169,5 +171,5 @@ _printf_(" configuring element and loads:\n"); - ConfigureObjectsx(elements, loads, nodes, materials,parameters); + ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters); _printf_(" process parameters:\n"); @@ -178,5 +180,5 @@ /*Use all the data created to create an femmodel: */ - femmodel=new FemModel(elements,nodes,constraints,loads,materials,parameters, + femmodel=new FemModel(elements,nodes,vertices, constraints,loads,materials,parameters, partition,tpartition,yg,Rmg,Gmn,nodesets,ys,ys0); Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 3416) +++ /issm/trunk/src/c/objects/Node.cpp (revision 3417) @@ -3,5 +3,5 @@ */ - +/*Include files: {{{1*/ #ifdef HAVE_CONFIG_H #include "config.h" @@ -11,5 +11,8 @@ #include "stdio.h" +#include "./Vertex.h" #include "./Node.h" +#include "./Hook.h" +#include "./DofIndexing.h" #include #include "../EnumDefinitions/EnumDefinitions.h" @@ -18,150 +21,162 @@ #include "../include/typedefs.h" #include "../include/macros.h" - -/*Object constructors and destructor*/ -/*FUNCTION Node constructor {{{1*/ +/*}}}*/ +/*Object constructors and destructors: {{{1*/ +/*FUNCTION Node default constructor {{{2*/ Node::Node(){ return; } /*}}}*/ -/*FUNCTION Node constructor {{{1*/ -Node::Node(int node_id,int node_partitionborder,int node_numdofs, double node_x[3],double node_sigma,int node_onbed,int node_onsurface,int node_upper_node_id,int node_onshelf,int node_onsheet){ - - int i; - - id=node_id; - partitionborder=node_partitionborder; - numberofdofs=node_numdofs; - x[0]=node_x[0]; - x[1]=node_x[1]; - x[2]=node_x[2]; - sigma=node_sigma; - onbed=node_onbed; - onsurface=node_onsurface; - onshelf=node_onshelf; - onsheet=node_onsheet; - - /*Initialize sets: */ - for(i=0;iid=node_id; return; } /*}}}*/ -/*FUNCTION Node destructor{{{1*/ +/*FUNCTION Node other constructor {{{2*/ +Node::Node(int node_id,DofIndexing* node_indexing, NodeProperties* node_properties, Hook* node_vertex, Hook* node_upper_node): + indexing(node_indexing), + properties(node_properties), + hvertex(node_vertex), + hupper_node(node_upper_node) { + + /*all the initialization has been done by the initializer, just fill in the id: */ + this->id=node_id; + +} +/*}}}*/ +/*FUNCTION Node constructor from iomodel{{{2*/ +Node::Node(int i, IoModel* iomodel){ //i is the node index + + int k; + + int numdofs; + int partitionborder; + int vertex_id; + + + /*id: */ + this->id=i+1; //matlab indexing + + /*indexing:*/ + DistributeNumDofs(&numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); //number of dofs per node + if(my_bordervertices[i])partitionborder=1; else partitionborder=0;//is this node on a partition border? + + this->indexing.Init(numdofs,partitionborder); + + /*properties: */ + this->properties.Init( + (int)iomodel->gridonbed[i], + (int)iomodel->gridonsurface[i], + (int)iomodel->gridoniceshelf[i], + (int)iomodel->gridonicesheet[i]); + + /*hooks: */ + vertex_id=this->id; //node and vertex have the same id, as we are running galerkin continuous, with same number of nodes and vertices. + + if (strcmp(iomodel->meshtype,"3d")==0){ + if (isnan(iomodel->uppernodes[i])){ + upper_node_id=this->id; //nodes on surface do not have upper nodes, only themselves. + } + else{ + upper_node_id=(int)iomodel->uppernodes[i]; + } + } + else{ + /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/ + upper_node_id=this->id; + } + + this->hvertex->Init(vertex_id,1); //node id is the same as the vertex id, continuous galerkin! + this->hupper_node->Init(upper_node_id,1); + + + /*set single point constraints.: */ + if (strcmp(iomodel->meshtype,"3d")==0){ + /*We have a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */ + if (iomodel->deadgrids[i]){ + for(k=1;k<=numdofs;k++){ + node->FreezeDof(k); + } + } + } + if (iomodel->gridonhutter[i]){ + for(k=1;k<=numdofs;k++){ + node->FreezeDof(k); + } + } +} +/*}}}*/ +/*FUNCTION Node destructor{{{2*/ Node::~Node(){ return; } /*}}}*/ - -/*Object management*/ -/*FUNCTION Node Marshall{{{1*/ -void Node::Marshall(char** pmarshalled_dataset){ +/*}}}*/ +/*Object management: {{{1*/ +/*FUNCTION Node Configure {{{2*/ +void Node::Configure(void* pnodes, void* pvertices){ + + int i; + + DataSet* nodesin=NULL; + DataSet* verticesin=NULL; + + /*Recover pointers :*/ + nodesin=(DataSet*)pnodes; + verticesin=(DataSet*)pvertices; + + /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective + * datasets, using internal ids and ofthis->indexing.f_sets hidden in hooks: */ + hvertex.configure(verticesin); + hupper_node.configure(nodesin); + +} +/*}}}*/ +/*FUNCTION Node copy {{{2*/ +Object* Node::copy() { + + return new Node(this->id,&this->indexing, &this->properties, &this->hvertex,&this->hupper_node); + +} + +/*}}}*/ +/*FUNCTION Node DeepEcho{{{2*/ +void Node::DeepEcho(void){ + + printf("Node:\n"); + printf(" id: %i\n",id); + indexing.DeepEcho(); + properties.DeepEcho(); + hvertex.DeepEcho(); + hupper_node.DeepEcho(); + +} +/*}}}*/ +/*FUNCTION Node Demarshall{{{2*/ +void Node::Demarshall(char** pmarshalled_dataset){ char* marshalled_dataset=NULL; - int enum_type=0; /*recover marshalled_dataset: */ marshalled_dataset=*pmarshalled_dataset; - /*get enum type of Node: */ - enum_type=NodeEnum(); - - /*marshall enum: */ - memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); - - /*marshall Node data: */ - memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id); - memcpy(marshalled_dataset,&partitionborder,sizeof(partitionborder));marshalled_dataset+=sizeof(partitionborder); - memcpy(marshalled_dataset,&clone,sizeof(clone));marshalled_dataset+=sizeof(clone); - memcpy(marshalled_dataset,&numberofdofs,sizeof(numberofdofs));marshalled_dataset+=sizeof(numberofdofs); - memcpy(marshalled_dataset,&x,sizeof(x));marshalled_dataset+=sizeof(x); - memcpy(marshalled_dataset,&sigma,sizeof(sigma));marshalled_dataset+=sizeof(sigma); - memcpy(marshalled_dataset,&onbed,sizeof(onbed));marshalled_dataset+=sizeof(onbed); - memcpy(marshalled_dataset,&onsurface,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface); - memcpy(marshalled_dataset,&onshelf,sizeof(onshelf));marshalled_dataset+=sizeof(onshelf); - memcpy(marshalled_dataset,&onsheet,sizeof(onsheet));marshalled_dataset+=sizeof(onsheet); - memcpy(marshalled_dataset,&doflist,sizeof(doflist));marshalled_dataset+=sizeof(doflist); - memcpy(marshalled_dataset,&doflist1,sizeof(doflist1));marshalled_dataset+=sizeof(doflist1); - memcpy(marshalled_dataset,&mset,sizeof(mset));marshalled_dataset+=sizeof(mset); - memcpy(marshalled_dataset,&nset,sizeof(nset));marshalled_dataset+=sizeof(nset); - memcpy(marshalled_dataset,&fset,sizeof(fset));marshalled_dataset+=sizeof(fset); - memcpy(marshalled_dataset,&sset,sizeof(sset));marshalled_dataset+=sizeof(sset); - memcpy(marshalled_dataset,&upper_node_id,sizeof(upper_node_id));marshalled_dataset+=sizeof(upper_node_id); - memcpy(marshalled_dataset,&upper_node,sizeof(upper_node));marshalled_dataset+=sizeof(upper_node); - memcpy(marshalled_dataset,&upper_node_offset,sizeof(upper_node_offset));marshalled_dataset+=sizeof(upper_node_offset); - - *pmarshalled_dataset=marshalled_dataset; - return; -} -/*}}}*/ -/*FUNCTION Node MarshallSize{{{1*/ -int Node::MarshallSize(){ - - return sizeof(id)+ - sizeof(partitionborder)+ - sizeof(clone)+ - sizeof(numberofdofs)+ - sizeof(x)+ - sizeof(sigma)+ - sizeof(onbed)+ - sizeof(onsurface)+ - sizeof(onshelf)+ - sizeof(onsheet)+ - sizeof(doflist)+ - sizeof(doflist1)+ - sizeof(mset)+ - sizeof(nset)+ - sizeof(fset)+ - sizeof(sset)+ - sizeof(upper_node_id)+ - sizeof(upper_node)+ - sizeof(upper_node_offset)+ - sizeof(int); //sizeof(int) for enum type -} -/*}}}*/ -/*FUNCTION Node Demarshall{{{1*/ -void Node::Demarshall(char** pmarshalled_dataset){ - - char* marshalled_dataset=NULL; - - /*recover marshalled_dataset: */ - marshalled_dataset=*pmarshalled_dataset; - /*this time, no need to get enum type, the pointer directly points to the beginning of the *object data (thanks to DataSet::Demarshall):*/ memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id); - memcpy(&partitionborder,marshalled_dataset,sizeof(partitionborder));marshalled_dataset+=sizeof(partitionborder); - memcpy(&clone,marshalled_dataset,sizeof(clone));marshalled_dataset+=sizeof(clone); - memcpy(&numberofdofs,marshalled_dataset,sizeof(numberofdofs));marshalled_dataset+=sizeof(numberofdofs); - memcpy(&x,marshalled_dataset,sizeof(x));marshalled_dataset+=sizeof(x); - memcpy(&sigma,marshalled_dataset,sizeof(sigma));marshalled_dataset+=sizeof(sigma); - memcpy(&onbed,marshalled_dataset,sizeof(onbed));marshalled_dataset+=sizeof(onbed); - memcpy(&onsurface,marshalled_dataset,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface); - memcpy(&onshelf,marshalled_dataset,sizeof(onshelf));marshalled_dataset+=sizeof(onshelf); - memcpy(&onsheet,marshalled_dataset,sizeof(onsheet));marshalled_dataset+=sizeof(onsheet); - memcpy(&doflist,marshalled_dataset,sizeof(doflist));marshalled_dataset+=sizeof(doflist); - memcpy(&doflist1,marshalled_dataset,sizeof(doflist1));marshalled_dataset+=sizeof(doflist1); - memcpy(&mset,marshalled_dataset,sizeof(mset));marshalled_dataset+=sizeof(mset); - memcpy(&nset,marshalled_dataset,sizeof(nset));marshalled_dataset+=sizeof(nset); - memcpy(&fset,marshalled_dataset,sizeof(fset));marshalled_dataset+=sizeof(fset); - memcpy(&sset,marshalled_dataset,sizeof(sset));marshalled_dataset+=sizeof(sset); - memcpy(&upper_node_id,marshalled_dataset,sizeof(upper_node_id));marshalled_dataset+=sizeof(upper_node_id); - memcpy(&upper_node,marshalled_dataset,sizeof(upper_node));marshalled_dataset+=sizeof(upper_node); - memcpy(&upper_node_offset,marshalled_dataset,sizeof(upper_node_offset));marshalled_dataset+=sizeof(upper_node_offset); - - /*upper node is not pointing to correct object anymore: */ - upper_node=NULL; - + + /*demarshall objects: */ + indexing.Demarshall(&marshalled_dataset); + properties.Demarshall(&marshalled_dataset); + hvertex.Demarshall(&marshalled_dataset); + hupper_node.Demarshall(&marshalled_dataset); + /*return: */ *pmarshalled_dataset=marshalled_dataset; @@ -169,15 +184,92 @@ } /*}}}*/ -/*FUNCTION Node GetId{{{1*/ +/*FUNCTION Node Echo{{{2*/ +void Node::Echo(void){ + + printf("Node:\n"); + printf(" id: %i\n",id); + indexing.Echo(); + properties.Echo(); + hvertex.Echo(); + hupper_node.Echo(); + +} +/*}}}*/ +/*FUNCTION Node Enum{{{2*/ +int Node::Enum(void){ + + return NodeEnum(); + +} +/*}}}*/ +/*FUNCTION Node GetId{{{2*/ int Node::GetId(void){ return id; } /*}}}*/ -/*FUNCTION Node GetName{{{1*/ +/*FUNCTION Node GetName{{{2*/ char* Node::GetName(void){ return "node"; } /*}}}*/ - -/*Object functions*/ -/*FUNCTION Node ApplyConstraints{{{1*/ +/*FUNCTION Node GetVertexDof {{{2*/ +int Node::GetVertexDof(void){ + + Vertex* vertex=NULL; + + vertex=(Vertex*)hvertex.delivers(); + return vertex->dof; +} +/*}}}*/ +/*FUNCTION Node Marshall{{{2*/ +void Node::Marshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int enum_type=0; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*get enum type of Node: */ + enum_type=NodeEnum(); + + /*marshall enum: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + + /*marshall Node data: */ + memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id); + + /*marshall objects: */ + indexing.Marshall(&marshalled_dataset); + properties.Marshall(&marshalled_dataset); + hvertex.Marshall(&marshalled_dataset); + hupper_node.Marshall(&marshalled_dataset); + + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ +/*FUNCTION Node MarshallSize{{{2*/ +int Node::MarshallSize(){ + + return sizeof(id)+ + +indexing.MarshallSize()+ + +properties.MarshallSize()+ + +hvertex.MarshallSize()+ + +hupper_node.MarshallSize()+ + sizeof(int); //sizeof(int) for enum type +} +/*}}}*/ +/*FUNCTION Node SetVertexDof {{{2*/ +void Node::SetVertexDof(int in_dof){ + + Vertex* vertex=NULL; + + vertex=(Vertex*)hvertex.delivers(); + vertex->dof=in_dof; + +} +/*}}}*/ +/*}}}*/ +/*Object numerics: {{{1*/ +/*FUNCTION Node ApplyConstraints{{{2*/ void Node::ApplyConstraint(Vec yg,int dof,double value){ @@ -192,7 +284,7 @@ * we are a clone!*/ - if(!clone){ - - index=doflist[dof-1]; //matlab indexing + if(!indexing.clone){ + + index=indexing.doflist[dof-1]; //matlab indexing VecSetValues(yg,1,&index,&value,INSERT_VALUES); @@ -202,23 +294,5 @@ } /*}}}*/ -/*FUNCTION Node Configure {{{1*/ -void Node::Configure(void* pnodes){ - - DataSet* nodes=NULL; - - /*Recover pointers :*/ - nodes=(DataSet*)pnodes; - - /*Link this node with its upper node: */ - ResolvePointers((Object**)&upper_node,&upper_node_id,&upper_node_offset,1,nodes); - -} -/*}}}*/ -/*FUNCTION Node copy {{{1*/ -Object* Node::copy() { - return new Node(*this); -} -/*}}}*/ -/*FUNCTION Node CreatePartition{{{1*/ +/*FUNCTION Node CreatePartition{{{2*/ void Node::CreatePartition(Vec partition){ @@ -227,5 +301,5 @@ idxm=(id-1); - value=(double)doflist1[0]; + value=(double)this->GetVertexDof(); ISSMASSERT(value>=0); @@ -235,5 +309,5 @@ } /*}}}*/ -/*FUNCTION Node CreateVecSets {{{1*/ +/*FUNCTION Node CreateVecSets {{{2*/ void Node::CreateVecSets(Vec pv_g,Vec pv_m,Vec pv_n,Vec pv_f,Vec pv_s){ @@ -243,69 +317,31 @@ int i; - for(i=0;iindexing.numberofdofs;i++){ /*g set: */ - VecSetValues(pv_g,1,&doflist[i],&gvalue,INSERT_VALUES); + VecSetValues(pv_g,1,&indexing.doflist[i],&gvalue,INSERT_VALUES); /*m set: */ - value=(double)mset[i]; - VecSetValues(pv_m,1,&doflist[i],&value,INSERT_VALUES); + value=(double)this->indexing.m_set[i]; + VecSetValues(pv_m,1,&indexing.doflist[i],&value,INSERT_VALUES); /*n set: */ - value=(double)nset[i]; - VecSetValues(pv_n,1,&doflist[i],&value,INSERT_VALUES); + value=(double)this->indexing.n_set[i]; + VecSetValues(pv_n,1,&indexing.doflist[i],&value,INSERT_VALUES); /*f set: */ - value=(double)fset[i]; - VecSetValues(pv_f,1,&doflist[i],&value,INSERT_VALUES); + value=(double)this->indexing.f_set[i]; + VecSetValues(pv_f,1,&indexing.doflist[i],&value,INSERT_VALUES); /*s set: */ - value=(double)sset[i]; - VecSetValues(pv_s,1,&doflist[i],&value,INSERT_VALUES); - - } - - -} -/*}}}*/ -/*FUNCTION Node DeepEcho{{{1*/ -void Node::DeepEcho(void){ - - int i; - - printf("Node:\n"); - printf(" id: %i\n",id); - printf(" partitionborder: %i\n",partitionborder); - printf(" clone: %i\n",clone); - printf(" numberofdofs: %i\n",numberofdofs); - printf(" x=[%g,%g,%g]\n",x[0],x[1],x[2]); - printf(" sigma=%g\n",sigma); - printf(" onbed: %i\n",onbed); - printf(" onsurface: %i\n",onsurface); - printf(" onshelf: %i\n",onshelf); - printf(" onsheet: %i\n",onsheet); - printf(" upper_node_id=%i\n",upper_node_id); - printf(" upper_node_offset=%i\n",upper_node_offset); - printf(" doflist:|"); - for(i=0;iMAXDOFSPERNODE)break; - printf("%i|",doflist[i]); - } - printf(" doflist1:|"); - printf("%i|",doflist1[0]); - printf("\n"); - - printf(" set membership: m,n,f,s sets \n"); - for(i=0;iMAXDOFSPERNODE)break; - printf(" dof %i: %i %i %i %i\n",i,mset[i],nset[i],fset[i],sset[i]); - } - printf("\n"); - if(upper_node)printf(" upper_node pointer: %p\n",upper_node); - - return; -} -/*}}}*/ -/*FUNCTION Node DistributeDofs{{{1*/ + value=(double)this->indexing.s_set[i]; + VecSetValues(pv_s,1,&indexing.doflist[i],&value,INSERT_VALUES); + + } + + +} +/*}}}*/ +/*FUNCTION Node DistributeDofs{{{2*/ void Node::DistributeDofs(int* pdofcount,int* pdofcount1){ @@ -318,5 +354,5 @@ dofcount1=*pdofcount1; - if(clone){ + if(indexing.clone){ /*This node is a clone! Don't distribute dofs, it will get them from another cpu!*/ return; @@ -324,10 +360,10 @@ /*This node should distribute dofs, go ahead: */ - for(i=0;iindexing.numberofdofs;i++){ + indexing.doflist[i]=dofcount+i; + } + dofcount+=this->indexing.numberofdofs; + + SetVertexDof(dofcount1); dofcount1+=1; @@ -339,16 +375,16 @@ } /*}}}*/ -/*FUNCTION Node DofInMSet{{{1*/ +/*FUNCTION Node DofInMSet{{{2*/ void Node::DofInMSet(int dof){ /*Put dof for this node into the m set (m set is for rigid body modes)*/ - mset[dof]=1; //m and n are mutually exclusive (m for rigid body modes) - nset[dof]=0; - fset[dof]=0; //n splits into f (for which we solve) and s (single point constraints) - sset[dof]=0; -} -/*}}}*/ -/*FUNCTION Node DofInSSet {{{1*/ + this->indexing.m_set[dof]=1; //m and n are mutually exclusive (m for rigid body modes) + this->indexing.n_set[dof]=0; + this->indexing.f_set[dof]=0; //n splits into f (for which we solve) and s (single point constraints) + this->indexing.s_set[dof]=0; +} +/*}}}*/ +/*FUNCTION Node DofInSSet {{{2*/ void Node::DofInSSet(int dof){ @@ -356,64 +392,19 @@ * to a fixed value during computations. */ - mset[dof]=0; //m and n are mutually exclusive (m for rigid body modes) - nset[dof]=1; - fset[dof]=0; //n splits into f (for which we solve) and s (single point constraints) - sset[dof]=1; -} -/*}}}*/ -/*FUNCTION Node DofIsInMSet{{{1*/ + this->indexing.m_set[dof]=0; //m and n are mutually exclusive (m for rigid body modes) + this->indexing.n_set[dof]=1; + this->indexing.f_set[dof]=0; //n splits into f (for which we solve) and s (single point constraints) + this->indexing.s_set[dof]=1; +} +/*}}}*/ +/*FUNCTION Node DofIsInMSet{{{2*/ int Node::DofIsInMSet(int dof){ - if (mset[dof])return 1; + if (this->indexing.m_set[dof])return 1; else return 0; } /*}}}*/ -/*FUNCTION Node Echo{{{1*/ -void Node::Echo(void){ - - int i; - - printf("Node:\n"); - printf(" id: %i\n",id); - printf(" partitionborder: %i\n",partitionborder); - printf(" clone: %i\n",clone); - printf(" numberofdofs: %i\n",numberofdofs); - printf(" x=[%g,%g,%g]\n",x[0],x[1],x[2]); - printf(" sigma=%g\n",sigma); - printf(" onbed: %i\n",onbed); - printf(" onsurface: %i\n",onsurface); - printf(" onshelf: %i\n",onshelf); - printf(" onsheet: %i\n",onsheet); - printf(" upper_node_id=%i\n",upper_node_id); - printf(" upper_node_offset=%i\n",upper_node_offset); - printf(" doflist:|"); - for(i=0;iMAXDOFSPERNODE)break; - printf("%i|",doflist[i]); - } - printf(" doflist1:|"); - printf("%i|",doflist1[0]); - printf("\n"); - - printf(" set membership: m,n,f,s sets \n"); - for(i=0;iMAXDOFSPERNODE)break; - printf(" dof %i: %i %i %i %i\n",i,mset[i],nset[i],fset[i],sset[i]); - } - printf("\n"); - if(upper_node)printf(" upper_node pointer: %p\n",upper_node); - - return; -} -/*}}}*/ -/*FUNCTION Node Enum{{{1*/ -int Node::Enum(void){ - - return NodeEnum(); - -} -/*}}}*/ -/*FUNCTION Node FieldAverageOntoVertices{{{1*/ +/*FUNCTION Node FieldAverageOntoVertices{{{2*/ void Node::FieldAverageOntoVertices(Vec field,double* field_serial,char* fieldname){ @@ -421,9 +412,9 @@ } /*}}}*/ -/*FUNCTION Node FieldDepthAverageAtBase{{{1*/ +/*FUNCTION Node FieldDepthAverageAtBase{{{2*/ void Node::FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname){ /* node data: */ - int doflist1; + int vertexdof; int dofx,dofy; int isnodeonsurface; @@ -436,7 +427,7 @@ /*Are we on the base, not on the surface, and not on a clone node?:*/ - if(onbed==1 & clone==0 &onsurface==0){ + if(properties.onbed==1 & indexing.clone==0 &properties.onsurface==0){ - doflist1=GetDofList1(); + vertexdof=this->GetVertexDof(); /*this node is on the bed. We are going to, follow the upper nodes until we reach the surface. At each upper node, @@ -456,6 +447,6 @@ /*get dofs for this base node velocity: we know there are two dofs in field_serial */ - dofx=2*doflist1; - dofy=2*doflist1+1; + dofx=2*vertexdof; + dofy=2*vertexdof+1; node=this; @@ -464,15 +455,15 @@ if (node->IsOnSurface())break; - doflist1=node->GetDofList1(); + vertexdof=node->GetVertexDof(); - velocity1[0]=field_serial[2*doflist1]; - velocity1[1]=field_serial[2*doflist1+1]; + velocity1[0]=field_serial[2*vertexdof]; + velocity1[1]=field_serial[2*vertexdof+1]; z1=node->GetZ(); upper_node=node->GetUpperNode(); - doflist1=upper_node->GetDofList1(); + vertexdof=upper_node->GetVertexDof(); - velocity2[0]=field_serial[2*doflist1]; - velocity2[1]=field_serial[2*doflist1+1]; + velocity2[0]=field_serial[2*vertexdof]; + velocity2[1]=field_serial[2*vertexdof+1]; z2=upper_node->GetZ(); @@ -507,5 +498,5 @@ /*get dofs for this base node velocity: we know there are two dofs in field_serial */ - dofx=doflist1; + dofx=vertexdof; node=this; @@ -514,13 +505,13 @@ if (node->IsOnSurface())break; - doflist1=node->GetDofList1(); + vertexdof=node->GetVertexDof(); - field1=field_serial[doflist1]; + field1=field_serial[vertexdof]; z1=node->GetZ(); upper_node=node->GetUpperNode(); - doflist1=upper_node->GetDofList1(); + vertexdof=upper_node->GetVertexDof(); - field2=field_serial[doflist1]; + field2=field_serial[vertexdof]; z2=upper_node->GetZ(); @@ -543,5 +534,5 @@ } /*}}}*/ -/*FUNCTION Node FieldExtrude {{{1*/ +/*FUNCTION Node FieldExtrude {{{2*/ void Node::FieldExtrude(Vec field,double* field_serial,char* field_name){ @@ -552,18 +543,18 @@ /*Is this node on bed? :*/ - if (onbed){ + if (properties.onbed){ if (strcmp(field_name,"velocity")==0){ /* node data: */ - int doflist1; + int vertexdof; int doflist[2]; double fieldnode[2]; - doflist1=GetDofList1(); + vertexdof=GetVertexDof(); /*get dofs for this base node velocity: we know there are two dofs in field_serial */ - doflist[0]=2*doflist1; - doflist[1]=2*doflist1+1; + doflist[0]=2*vertexdof; + doflist[1]=2*vertexdof+1; //initilaize node @@ -578,7 +569,7 @@ for(;;){ - doflist1=node->GetDofList1(); - doflist[0]=2*doflist1; - doflist[1]=2*doflist1+1; + vertexdof=node->GetVertexDof(); + doflist[0]=2*vertexdof; + doflist[1]=2*vertexdof+1; VecSetValues(field,1,&doflist[0],&fieldnode[0],INSERT_VALUES); VecSetValues(field,1,&doflist[1],&fieldnode[1],INSERT_VALUES); @@ -597,5 +588,5 @@ //initilaize node and get dof1 node=this; - dof1=node->GetDofList1(); + dof1=node->GetVertexDof(); /*get field for this base node: */ @@ -606,5 +597,5 @@ for(;;){ - dof1=node->GetDofList1(); + dof1=node->GetVertexDof(); VecSetValues(field,1,&dof1,&fieldel,INSERT_VALUES); @@ -655,5 +646,5 @@ } /*}}}*/ -/*FUNCTION Node FreezeDof{{{1*/ +/*FUNCTION Node FreezeDof{{{2*/ void Node::FreezeDof(int dof){ @@ -662,79 +653,96 @@ } /*}}}*/ -/*FUNCTION Node GetDof {{{1*/ +/*FUNCTION Node GetDof {{{2*/ int Node::GetDof(int dofindex){ - return doflist[dofindex]; - -} -/*}}}*/ -/*FUNCTION Node GetDofList{{{1*/ + return indexing.doflist[dofindex]; + +} +/*}}}*/ +/*FUNCTION Node GetDofList{{{2*/ void Node::GetDofList(int* outdoflist,int* pnumberofdofspernode){ int i; - for(i=0;iindexing.numberofdofs;i++){ + outdoflist[i]=indexing.doflist[i]; } /*Assign output pointers:*/ - *pnumberofdofspernode=numberofdofs; -} -/*}}}*/ -/*FUNCTION Node GetDOfList1{{{1*/ -int Node::GetDofList1(void){ - return doflist1[0]; -} -/*}}}*/ -/*FUNCTION Node GetNumberOfDofs{{{1*/ + *pnumberofdofspernode=this->indexing.numberofdofs; +} +/*}}}*/ +/*FUNCTION Node GetNumberOfDofs{{{2*/ int Node::GetNumberOfDofs(){ - return numberofdofs; - -} -/*}}}*/ -/*FUNCTION Node GetSigma {{{1*/ -double Node::GetSigma(){return sigma;} -/*}}}*/ -/*FUNCTION Node GetUpperNode {{{1*/ + return this->indexing.numberofdofs; + +} +/*}}}*/ +/*FUNCTION Node GetSigma {{{2*/ +double Node::GetSigma(){ + Vertex* vertex=NULL; + + vertex=(Vertex*)hvertex.delivers(); + return vertex->sigma; +} +/*}}}*/ +/*FUNCTION Node GetUpperNode {{{2*/ Node* Node::GetUpperNode(){ + Node* upper_node=NULL; + upper_node=(Node*)hupper_node.delivers(); return upper_node; } /*}}}*/ -/*FUNCTION Node GetX {{{1*/ -double Node::GetX(){return x[0];} -/*}}}*/ -/*FUNCTION Node GetY {{{1*/ -double Node::GetY(){return x[1];} -/*}}}*/ -/*FUNCTION Node GetZ {{{1*/ -double Node::GetZ(){return x[2];} -/*}}}*/ -/*FUNCTION Node IsClone {{{1*/ +/*FUNCTION Node GetX {{{2*/ +double Node::GetX(){ + Vertex* vertex=NULL; + + vertex=(Vertex*)hvertex.delivers(); + return vertex->x; +} +/*}}}*/ +/*FUNCTION Node GetY {{{2*/ +double Node::GetY(){ + Vertex* vertex=NULL; + + vertex=(Vertex*)hvertex.delivers(); + return vertex->y; +} +/*}}}*/ +/*FUNCTION Node GetZ {{{2*/ +double Node::GetZ(){ + Vertex* vertex=NULL; + + vertex=(Vertex*)hvertex.delivers(); + return vertex->z; +} +/*}}}*/ +/*FUNCTION Node IsClone {{{2*/ int Node::IsClone(){ - return clone; - -} -/*}}}*/ -/*FUNCTION Node IsOnBed {{{1*/ + return indexing.clone; + +} +/*}}}*/ +/*FUNCTION Node IsOnBed {{{2*/ int Node::IsOnBed(){ - return onbed; -} -/*}}}*/ -/*FUNCTION Node IsOnSheet {{{1*/ + return properties.onbed; +} +/*}}}*/ +/*FUNCTION Node IsOnSheet {{{2*/ int Node::IsOnSheet(){ - return onsheet; + return properties.onsheet; } /*}}}*/ -/*FUNCTION Node IsOnShelf {{{1*/ +/*FUNCTION Node IsOnShelf {{{2*/ int Node::IsOnShelf(){ - return onshelf; -} -/*}}}*/ -/*FUNCTION Node IsOnSurface {{{1*/ + return properties.onshelf; +} +/*}}}*/ +/*FUNCTION Node IsOnSurface {{{2*/ int Node::IsOnSurface(){ - return onsurface; -} -/*}}}*/ -/*FUNCTION Node MyRank{{{1*/ + return properties.onsurface; +} +/*}}}*/ +/*FUNCTION Node MyRank{{{2*/ int Node::MyRank(void){ extern int my_rank; @@ -743,5 +751,5 @@ } /*}}}*/ -/*FUNCTION Node SetClone {{{1*/ +/*FUNCTION Node SetClone {{{2*/ void Node::SetClone(int* minranks){ @@ -749,15 +757,15 @@ if (minranks[id-1]==my_rank){ - clone=0; + indexing.clone=0; } else{ /*!there is a cpu with lower rank that has the same node, therefore, I am a clone*/ - clone=1; - } - -} -/*}}}*/ -/*FUNCTION Node ShowBorderDofs{{{1*/ + indexing.clone=1; + } + +} +/*}}}*/ +/*FUNCTION Node ShowBorderDofs{{{2*/ void Node::ShowBorderDofs(int* borderdofs,int* borderdofs1){ @@ -766,20 +774,20 @@ /*Is this node on the partition border? */ - if(!partitionborder)return; + if(!indexing.partitionborder)return; /*Are we the cpu that created this node's dof list? */ - if(clone)return; + if(indexing.clone)return; /*Ok, we are on the partition border, and we did create the * dofs for this node, plug the doflist into borderdofs: */ - for(j=0;jindexing.numberofdofs;j++){ + *(borderdofs+this->indexing.numberofdofs*(id-1)+j)=indexing.doflist[j]; + } + *(borderdofs1+(id-1)+0)=this->GetVertexDof(); return; } /*}}}*/ -/*FUNCTION Node UpdateBorderDofs{{{1*/ +/*FUNCTION Node UpdateBorderDofs{{{2*/ void Node::UpdateBorderDofs(int* allborderdofs,int* allborderdofs1){ @@ -788,20 +796,20 @@ /*Is this node on the partition border? */ - if(!partitionborder)return; + if(!indexing.partitionborder)return; /*Are we the cpu that created this node's dof list? */ - if(clone==0)return; + if(indexing.clone==0)return; /*Ok, we are on the partition border, but we did not create * the dofs for this node. Therefore, our doflist is garbage right * now. Go pick it up in the allborderdofs: */ - for(j=0;jindexing.numberofdofs;j++){ + indexing.doflist[j]=*(allborderdofs+this->indexing.numberofdofs*(id-1)+j); + } + this->SetVertexDof(*(allborderdofs1+(id-1)+0)); return; } /*}}}*/ -/*FUNCTION Node UpdateDofs{{{1*/ +/*FUNCTION Node UpdateDofs{{{2*/ void Node::UpdateDofs(int dofcount,int dofcount1){ @@ -809,5 +817,5 @@ extern int my_rank; - if(clone){ + if(indexing.clone){ /*This node is a clone, don't update the dofs!: */ return; @@ -815,40 +823,18 @@ /*This node should update the dofs, go ahead: */ - for(i=0;iindexing.numberofdofs;i++){ + indexing.doflist[i]+=dofcount; + } + this->SetVertexDof(this->GetVertexDof()+dofcount1); return; } /*}}}*/ -/*FUNCTION Node UpdateFromInputs {{{1*/ +/*FUNCTION Node UpdateFromInputs {{{2*/ void Node::UpdateFromInputs(void* vinputs){ - - ParameterInputs* inputs=NULL; - Node* node=this; - int dof[1]={0}; - - /*Recover parameter inputs: */ - inputs=(ParameterInputs*)vinputs; - - /*Update internal data if inputs holds new values: */ - inputs->Recover("x",&x[0],1,dof,1,(void**)&node); - inputs->Recover("y",&x[1],1,dof,1,(void**)&node); - inputs->Recover("z",&x[2],1,dof,1,(void**)&node); - - -} -/*}}}*/ -/*FUNCTION Node UpdateNodePosition {{{1*/ -void Node::UpdateNodePosition(double* thickness,double* bed){ - - int dof; - - dof=this->GetDofList1(); - - /*sigma remains constant. z=bed+sigma*thickness*/ - this->x[2]=bed[dof]+sigma*thickness[dof]; - -} -/*}}}*/ + + ISSMERROR("not used yet!"); + +} +/*}}}*/ +/*}}}*/ Index: /issm/trunk/src/c/objects/Node.h =================================================================== --- /issm/trunk/src/c/objects/Node.h (revision 3416) +++ /issm/trunk/src/c/objects/Node.h (revision 3417) @@ -6,55 +6,53 @@ #define _NODE_H_ +/*indefinitions: */ +class Object; +class Vertex; +class Hook; + #include "./Object.h" +#include "../DataSet/DataSet.h" +#include "./Hook.h" +#include "./Vertex.h" +#include "./DofIndexing.h" +#include "./NodeProperties.h" #include "../toolkits/toolkits.h" - -#define MAXDOFSPERNODE 4 -#define PSSTRINGLENGTH 20 class Node: public Object{ private: - /*raw data coming from the model processor: */ - int id; /*! id, to track it*/ - int partitionborder; /*! during parallel partitioning, does this grid belong to a partition border?*/ - int numberofdofs; //real number of dofs, between 0 and MAXDOFSPERNODE - int clone; /*!this nodes is one the partition border, and is cloned*/ - double x[3]; /*! coordinates*/ - double sigma; /*! sigma = (z-bed)/thickness*/ + + int id; + + DofIndexing indexing; + NodeProperties properties; + Hook hvertex; + Hook hupper_node; + - int onbed; /*! for 3d, on bedrock*/ - int onsurface; /*! for 3d, on surface*/ - int onshelf; - int onsheet; - - /*for dof constraining: */ - int mset[MAXDOFSPERNODE]; - int nset[MAXDOFSPERNODE]; - int fset[MAXDOFSPERNODE]; - int sset[MAXDOFSPERNODE]; - - int upper_node_id; - Node* upper_node; - int upper_node_offset; - - /*data that is post processed : */ - int doflist[MAXDOFSPERNODE]; //dof list on which we solve - int doflist1[1]; //1d dof list to recover input parameters - public: + /*FUNCTION constructors, destructors {{{1*/ Node(); - Node(int node_id,int node_partitionborder,int node_numdofs, double node_x[3],double sigma,int node_onbed,int node_onsurface,int upper_node_id,int onshelf,int onsheet); + Node(int id,int vertex_id, int upper_node_id, int partitionborder,int numberofdofs, NodeProperties* node_properties); + Node(int id,DofIndexing* indexing, NodeProperties* properties, Hook* vertex, Hook* upper_node); + Node(int i, IoModel* iomodel); ~Node(); - + /*}}}*/ + /*FUNCTION object management {{{1*/ + void Configure(void* pnodes, void* pvertices); + void DeepEcho(); + void Demarshall(char** pmarshalled_dataset); void Echo(); - void DeepEcho(); + int Enum(); + int GetId(void); + char* GetName(); + int GetVertexDof(void); void Marshall(char** pmarshalled_dataset); int MarshallSize(); - char* GetName(); - void Demarshall(char** pmarshalled_dataset); - int Enum(); - int GetId(void); int MyRank(void); + void SetVertexDof(int in_dof); + /*}}}*/ + /*FUNCTION numerical routines {{{1*/ void DistributeDofs(int* pdofcount,int* pdofcount1); void UpdateDofs(int dofcount,int dofcount1); @@ -79,5 +77,4 @@ Object* copy(); void UpdateFromInputs(void* inputs); - void Configure(void* pnodes); Node* GetUpperNode(); int IsOnBed(); @@ -90,6 +87,6 @@ void FieldExtrude(Vec field,double* field_serial,char* field_name); void UpdateNodePosition(double* thickness,double* bed); + /*}}}*/ }; #endif /* _NODE_H_ */ - Index: /issm/trunk/src/c/objects/NodeProperties.cpp =================================================================== --- /issm/trunk/src/c/objects/NodeProperties.cpp (revision 3417) +++ /issm/trunk/src/c/objects/NodeProperties.cpp (revision 3417) @@ -0,0 +1,132 @@ +/*!\file NodeProperties.c + * \brief: implementation of the NodeProperties object + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "stdio.h" +#include "./NodeProperties.h" +#include +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../include/typedefs.h" +#include "../include/macros.h" + +/*Object constructors and destructor*/ +/*FUNCTION NodeProperties default constructor {{{1*/ +NodeProperties::NodeProperties(){ + + /*initialize to UNDEF every value: */ + this->onbed=UNDEF; + this->onsurface=UNDEF; + this->onshelf=UNDEF; + this->onsheet=UNDEF; +} +/*}}}*/ +/*FUNCTION NodeProperties constructor {{{1*/ +NodeProperties::NodeProperties(int nodeproperties_onbed, int nodeproperties_onsurface, int nodeproperties_onshelf, int nodeproperties_onsheet){ + this->Init(nodeproperties_onbed, nodeproperties_onsurface, nodeproperties_onshelf, nodeproperties_onsheet); + +} +/*}}}*/ +/*FUNCTION NodeProperties init: used by constructor {{{1*/ +NodeProperties::Init(int nodeproperties_onbed, int nodeproperties_onsurface, int nodeproperties_onshelf, int nodeproperties_onsheet){ + + this->onbed=nodeproperties_onbed; + this->onsurface=nodeproperties_onsurface; + this->onshelf=nodeproperties_onshelf; + this->onsheet=nodeproperties_onsheet; + +} +/*}}}*/ +/*FUNCTION NodeProperties copy constructor{{{1*/ +NodeProperties::NodeProperties(NodeProperties* prop){ //copy constructor + + this->onbed=prop->onbed; + this->onsurface=prop->onsurface; + this->onshelf=prop->onshelf; + this->onsheet=prop->onsheet; + +} +/*}}}*/ +/*FUNCTION NodeProperties destructor {{{1*/ +NodeProperties::~NodeProperties(){ //destructor +} +/*}}}*/ +/*FUNCTION NodeProperties Marshall{{{1*/ +void NodeProperties::Marshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int enum_type=0; + int nill=0; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*get enum type of NodeProperties: */ + enum_type=NodePropertiesEnum(); + + /*marshall enum: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + + memcpy(marshalled_dataset,&onbed,sizeof(onbed));marshalled_dataset+=sizeof(onbed); + memcpy(marshalled_dataset,&onsurface,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface); + memcpy(marshalled_dataset,&onshelf,sizeof(onshelf));marshalled_dataset+=sizeof(onshelf); + memcpy(marshalled_dataset,&onsheet,sizeof(onsheet));marshalled_dataset+=sizeof(onsheet); + + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ +/*FUNCTION NodeProperties MarshallSize{{{1*/ +int NodeProperties::MarshallSize(){ + + return sizeof(onbed)+ + sizeof(onsurface)+ + sizeof(onshelf)+ + sizeof(onsheet)+ + sizeof(int); //sizeof(int) for enum type + +} +/*}}}*/ +/*FUNCTION NodeProperties Demarshall{{{1*/ +void NodeProperties::Demarshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int i; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*this time, no need to get enum type, the pointer directly points to the beginning of the + *object data (thanks to DataSet::Demarshall):*/ + + memcpy(&onbed,marshalled_dataset,sizeof(onbed));marshalled_dataset+=sizeof(onbed); + memcpy(&onsurface,marshalled_dataset,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface); + memcpy(&onshelf,marshalled_dataset,sizeof(onshelf));marshalled_dataset+=sizeof(onshelf); + memcpy(&onsheet,marshalled_dataset,sizeof(onsheet));marshalled_dataset+=sizeof(onsheet); + + /*return: */ + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*FUNCTION NodeProperties Echo{{{1*/ +void NodeProperties::Echo(void){ + + printf("NodeProperties echo: \n"); + printf(" onbed: %i\n",onbed); + printf(" onsurface: %i\n",onsurface); + printf(" onshelf: %i\n",onshelf); + printf(" onsheet: %i\n",onsheet); +} +/*}}}*/ +/*FUNCTION NodeProperties DeepEcho{{{1*/ +void NodeProperties::DeepEcho(void){ + this->Echo(); +} +/*}}}*/ Index: /issm/trunk/src/c/objects/NodeProperties.h =================================================================== --- /issm/trunk/src/c/objects/NodeProperties.h (revision 3417) +++ /issm/trunk/src/c/objects/NodeProperties.h (revision 3417) @@ -0,0 +1,33 @@ +/*!\file: NodeProperties.h + * \brief prototype for NodeProperties.h + */ + +#ifndef _NODEPROPERTIES_H_ +#define _NODEPROPERTIES_H_ + +class NodeProperties{ + + public: + + int onbed; /*! for 3d, on bedrock*/ + int onsurface; /*! for 3d, on surface*/ + int onshelf; + int onsheet; + + NodeProperties(); + NodeProperties(int onbed, int onsurface, int onshelf, int onsheet); + Init(int onbed, int onsurface, int onshelf, int onsheet); + NodeProperties(NodeProperties* properties); + ~NodeProperties(); + + void Marshall(char** pmarshalled_dataset); + int MarshallSize(); + void Demarshall(char** pmarshalled_dataset); + void copy(NodeProperties* properties); + void Echo(void); + void DeepEcho(void); + +}; +#endif //ifndef _NODEPROPERTIES_H_ + + Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 3416) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 3417) @@ -51,4 +51,64 @@ return; +} +/*}}}*/ +/*FUNCTION Penta other constructor {{{1*/ +Penta::Penta(int i, IoModel* iomodel){ //i is the element index + + int j; + int offset; + int penta_node_ids[6]; + double penta_h[6]; + double penta_s[6]; + double penta_b[6]; + double penta_k[6]; + double penta_melting[6]; + double penta_accumulation[6]; + + /*id: */ + this->id=i+1; + + /*hooks: */ + for(j=0;j<6;j++){ //go recover node ids, needed to initialize the node hook. + penta_node_ids[j]=(int)*(iomodel->elements+6*i+j); //ids for vertices are in the elements array from Matlab + } + penta_matice_id=i+1; //refers to the corresponding ice material object + penta_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object + penta_numpar_id=1; //refers to numerical parameters object + + this->hnodes->Init(penta_node_ids,6); + this->hmatice->Init(penta_matice_id,1); + this->hmatpar->Init(penta_matpar_id,1); + this->hnumpar->Init(penta_numpar_id,1); + + /*properties: */ + for(j=0;j<6;j++){ + + offset=(int)*(iomodel->elements+6*i+j)-1; //get index of this vertex, in C numbering. + penta_h[j]= *(iomodel->thickness+offset); + penta_s[j]= *(iomodel->surface+offset); + penta_b[j]=*(iomodel->bed+offset); + penta_k[j]=*(iomodel->drag+offset); + penta_melting[j]=*(iomodel->melting+offset); + penta_accumulation[j]=*(iomodel->accumulation+offset); + } + penta_friction_type=(int)iomodel->drag_type; + penta_p=iomodel->p[i]; + penta_q=iomodel->q[i]; + penta_shelf=(int)*(iomodel->elementoniceshelf+i); + penta_onwater=(bool)*(iomodel->elementonwater+i); + penta_onsurface=(int)*(iomodel->elementonsurface+i); + penta_onwater=(bool)*(iomodel->elementonwater+i); + + if (*(iomodel->elements_type+2*i+0)==macayealformulationenum()){ //elements of type 3 are macayeal type penta. we collapse the formulation on their base. + penta_collapse=1; + } + else{ + penta_collapse=0; + } + + this->properties->Init(6,penta_h, penta_s, penta_b, penta_k, penta_melting, penta_accumulation, NULL, penta_friction_type, penta_p, penta_q, penta_shelf, penta_onbed, penta_onwater, penta_onsurface, penta_collapse, UNDEF); + + } /*}}}*/ Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 3416) +++ /issm/trunk/src/c/objects/Penta.h (revision 3417) @@ -11,15 +11,15 @@ #include "./Numpar.h" #include "./Matice.h" -#include "./Node.h" #include "./Tria.h" #include "./Hook.h" #include "./ElementProperties.h" #include "./ParameterInputs.h" +#include "./Node.h" -class Numpar; +class Object; class Node; -class Matice; -class Matpar; -class ElementProperties; +class Hook; +class ElementProperties; +class DataSet; class Penta: public Element{ @@ -41,4 +41,5 @@ Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id, int penta_numpar_id, ElementProperties* properties); Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Hook* penta_hnumpar, ElementProperties* penta_properties); + Penta(int i, IoModel* iomodel); ~Penta(); /*}}}*/ Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 3416) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 3417) @@ -59,4 +59,56 @@ return; +} +/*}}}*/ +/*FUNCTION Tria constructor from iomodel{{{1*/ +Tria::Tria(int i, IoModel* iomodel){ //i is the element index + + int j; + int offset; + 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]; + + /*id: */ + this->id=i+1; + + /*hooks: */ + for(j=0;j<3;j++){ //go recover node ids, needed to initialize the node hook. + tria_node_ids[j]=(int)*(iomodel->elements+3*i+j); //ids for vertices are in the elements array from Matlab + } + tria_matice_id=i+1; //refers to the corresponding ice material object + tria_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object + tria_numpar_id=1; //refers to numerical parameters object + + this->hnodes->Init(tria_node_ids,3); + this->hmatice->Init(tria_matice_id,1); + this->hmatpar->Init(tria_matpar_id,1); + this->hnumpar->Init(tria_numpar_id,1); + + /*properties: */ + for(j=0;j<3;j++){ + + offset=(int)*(iomodel->elements+3*i+j)-1; //get index of this vertex, in C numbering. + tria_h[j]= *(iomodel->thickness+offset); + tria_s[j]= *(iomodel->surface+offset); + tria_b[j]=*(iomodel->bed+offset); + tria_k[j]=*(iomodel->drag+offset); + tria_melting[j]=*(iomodel->melting+offset); + tria_accumulation[j]=*(iomodel->accumulation+offset); + } + tria_friction_type=(int)iomodel->drag_type; + tria_p=iomodel->p[i]; + tria_q=iomodel->q[i]; + tria_shelf=(int)*(iomodel->elementoniceshelf+i); + tria_onwater=(bool)*(iomodel->elementonwater+i); + + this->properties->Init(3,tria_h, tria_s, tria_b, tria_k, tria_melting, tria_accumulation, NULL, + tria_friction_type, tria_p, tria_q, tria_shelf, UNDEF,tria_onwater, UNDEF,UNDEF,UNDEF); + + } /*}}}*/ Index: /issm/trunk/src/c/objects/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Tria.h (revision 3416) +++ /issm/trunk/src/c/objects/Tria.h (revision 3417) @@ -7,16 +7,17 @@ +#include "./Object.h" #include "./Element.h" +#include "./Hook.h" +#include "./Node.h" +#include "./ElementProperties.h" #include "../DataSet/DataSet.h" -#include "./Object.h" -#include "./Node.h" -#include "./Matice.h" -#include "./Matpar.h" -#include "./Numpar.h" -#include "./ParameterInputs.h" -#include "./ElementProperties.h" -#include "./Hook.h" class Object; +class Element; +class Hook; +class ElementProperties; +class DataSet; +class Node; class Tria: public Element{ @@ -39,4 +40,5 @@ Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id, int tria_numpar_id, ElementProperties* tria_properties); Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Hook* tria_hnumpar, ElementProperties* tria_properties); + Tria(int i, IoModel* iomodel); ~Tria(); /*}}}*/ Index: /issm/trunk/src/c/objects/Vertex.cpp =================================================================== --- /issm/trunk/src/c/objects/Vertex.cpp (revision 3416) +++ /issm/trunk/src/c/objects/Vertex.cpp (revision 3417) @@ -3,4 +3,5 @@ */ +/*Include files: {{{1*/ #ifdef HAVE_CONFIG_H #include "config.h" @@ -16,4 +17,5 @@ #include "../include/typedefs.h" #include "../include/macros.h" +/*}}}*/ /*Object constructors and destructor*/ @@ -24,5 +26,10 @@ /*}}}*/ /*FUNCTION Vertex constructor {{{1*/ -Vertex::Vertex(int tria_id, double tria_x, double tria_y, double tria_z){ +Vertex::Vertex(int tria_id, double tria_x, double tria_y, double tria_z, double tria_sigma, int partitionborder){ + this->Init(tria_id, tria_x, tria_y, tria_z, tria_sigma, partitionborder); +} +/*}}}*/ +/*FUNCTION Vertex init, used by constructor {{{1*/ +Vertex::Init(int tria_id, double tria_x, double tria_y, double tria_z, double tria_sigma, int partitionborder){ /*all the initialization has been done by the initializer, just fill in the id: */ @@ -31,6 +38,21 @@ this->y=tria_y; this->z=tria_z; - - return; + this->sigma=tria_sigma; + this->dof=UNDEF; + + return; +} +/*}}}*/ +/*FUNCTION Vertex constructor from iomodel{{{1*/ +Vertex::Vertex(int i, IoModel* iomodel){ + + int partitionborder; + + /*is this vertex on a partition border? */ + if(my_bordervertices[i])partitionborder=1; + else partitionborder=0; + + this->Init(i+1, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]),partitionborder); + } /*}}}*/ @@ -71,4 +93,8 @@ memcpy(&y,marshalled_dataset,sizeof(y));marshalled_dataset+=sizeof(y); memcpy(&z,marshalled_dataset,sizeof(z));marshalled_dataset+=sizeof(z); + memcpy(&sigma,marshalled_dataset,sizeof(sigma));marshalled_dataset+=sizeof(sigma); + memcpy(&dof,marshalled_dataset,sizeof(dof));marshalled_dataset+=sizeof(dof); + memcpy(&partitionborder,marshalled_dataset,sizeof(partitionborder));marshalled_dataset+=sizeof(partitionborder); + memcpy(&clone,marshalled_dataset,sizeof(clone));marshalled_dataset+=sizeof(clone); /*return: */ @@ -86,4 +112,8 @@ printf(" y: %g\n",y); printf(" z: %g\n",z); + printf(" sigma: %g\n",sigma); + printf(" dof: %g\n",dof); + printf(" partitionborder: %g\n",partitionborder); + printf(" clone: %g\n",clone); return; @@ -125,4 +155,8 @@ memcpy(marshalled_dataset,&y,sizeof(y));marshalled_dataset+=sizeof(y); memcpy(marshalled_dataset,&z,sizeof(z));marshalled_dataset+=sizeof(z); + memcpy(marshalled_dataset,&sigma,sizeof(sigma));marshalled_dataset+=sizeof(sigma); + memcpy(marshalled_dataset,&dof,sizeof(dof));marshalled_dataset+=sizeof(dof); + memcpy(marshalled_dataset,&partitionborder,sizeof(partitionborder));marshalled_dataset+=sizeof(partitionborder); + memcpy(marshalled_dataset,&clone,sizeof(clone));marshalled_dataset+=sizeof(clone); *pmarshalled_dataset=marshalled_dataset; @@ -137,4 +171,8 @@ sizeof(y)+ sizeof(z)+ + sizeof(sigma)+ + sizeof(dof)+ + sizeof(partitionborder)+ + sizeof(clone)+ +sizeof(int); //sizeof(int) for enum type } @@ -155,7 +193,32 @@ /*FUNCTION UpdateFromInputs {{{1*/ void Vertex::UpdateFromInputs(void* vinputs){ - + + ParameterInputs* inputs=NULL; + Vertex* vertex=NULL; + int dof[1]={0}; + double coord[3]; + + vertex=this; + + coord[0]=this->x; coord[1]=this->y; coord[2]=this->z; + + /*Recover parameter inputs: */ + inputs=(ParameterInputs*)vinputs; + + /*Update internal data if inputs holds new values: */ + inputs->Recover("x",&coord[0],1,dof,1,(void**)&vertex); + inputs->Recover("y",&coord[1],1,dof,1,(void**)&vertex); + inputs->Recover("z",&coord[2],1,dof,1,(void**)&vertex); + ISSMERROR("not supported yet!"); } /*}}}*/ +/*FUNCTION UpdateVertexPosition {{{1*/ +void Vertex::UpdatePosition(double* thickness,double* bed){ + + /*sigma remains constant. z=bed+sigma*thickness*/ + this->z=bed[dof]+sigma*thickness[this->dof]; + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Vertex.h =================================================================== --- /issm/trunk/src/c/objects/Vertex.h (revision 3416) +++ /issm/trunk/src/c/objects/Vertex.h (revision 3417) @@ -10,5 +10,5 @@ class Vertex: public Object{ - private: + public: int id; @@ -16,10 +16,16 @@ double y; double z; + double sigma; //sigma coordinate: (z-bed)/thickness - public: + /*dof management: */ + int partitionborder; + int clone; + int dof; //dof to recover values in a vertex indexed vector /*FUNCTION constructors, destructors {{{1*/ Vertex(); - Vertex(int id, double x, double y, double z); + Vertex(int id, double x, double y, double z, double sigma,int partitionborder); + Init(int id, double x, double y, double z, double sigma,int partitionborder); + Vertex(int i, IoModel* iomodel); ~Vertex(); /*}}}*/ @@ -37,4 +43,7 @@ void UpdateFromDakota(void* vinputs); void UpdateFromInputs(void* vinputs); + void UpdatePosition(double* thickness,double* bed); + + /*}}}*/ Index: /issm/trunk/src/m/enum/BeamEnum.m =================================================================== --- /issm/trunk/src/m/enum/BeamEnum.m (revision 3416) +++ /issm/trunk/src/m/enum/BeamEnum.m (revision 3417) @@ -7,3 +7,3 @@ % macro=BeamEnum() -macro=415; +macro=416; Index: /issm/trunk/src/m/enum/DofIndexingEnum.m =================================================================== --- /issm/trunk/src/m/enum/DofIndexingEnum.m (revision 3417) +++ /issm/trunk/src/m/enum/DofIndexingEnum.m (revision 3417) @@ -0,0 +1,9 @@ +function macro=DofIndexingEnum() +%DOFINDEXINGENUM - Enum of DofIndexing +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=DofIndexingEnum() + +macro=417; Index: /issm/trunk/src/m/enum/NodeEnum.m =================================================================== --- /issm/trunk/src/m/enum/NodeEnum.m (revision 3416) +++ /issm/trunk/src/m/enum/NodeEnum.m (revision 3417) @@ -7,3 +7,3 @@ % macro=NodeEnum() -macro=420; +macro=421; Index: /issm/trunk/src/m/enum/NodePropertiesEnum.m =================================================================== --- /issm/trunk/src/m/enum/NodePropertiesEnum.m (revision 3417) +++ /issm/trunk/src/m/enum/NodePropertiesEnum.m (revision 3417) @@ -0,0 +1,9 @@ +function macro=NodePropertiesEnum() +%NODEPROPERTIESENUM - Enum of NodeProperties +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=NodePropertiesEnum() + +macro=413; Index: /issm/trunk/src/m/enum/PentaEnum.m =================================================================== --- /issm/trunk/src/m/enum/PentaEnum.m (revision 3416) +++ /issm/trunk/src/m/enum/PentaEnum.m (revision 3417) @@ -7,3 +7,3 @@ % macro=PentaEnum() -macro=413; +macro=414; Index: /issm/trunk/src/m/enum/SingEnum.m =================================================================== --- /issm/trunk/src/m/enum/SingEnum.m (revision 3416) +++ /issm/trunk/src/m/enum/SingEnum.m (revision 3417) @@ -7,3 +7,3 @@ % macro=SingEnum() -macro=414; +macro=415; Index: /issm/trunk/src/m/enum/VertexEnum.m =================================================================== --- /issm/trunk/src/m/enum/VertexEnum.m (revision 3417) +++ /issm/trunk/src/m/enum/VertexEnum.m (revision 3417) @@ -0,0 +1,9 @@ +function macro=VertexEnum() +%VERTEXENUM - Enum of Vertex +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=VertexEnum() + +macro=420; Index: /issm/trunk/src/m/enum/VerticesEnum.m =================================================================== --- /issm/trunk/src/m/enum/VerticesEnum.m (revision 3417) +++ /issm/trunk/src/m/enum/VerticesEnum.m (revision 3417) @@ -0,0 +1,9 @@ +function macro=VerticesEnum() +%VERTICESENUM - Enum of Vertices +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=VerticesEnum() + +macro=108; Index: /issm/trunk/src/m/solutions/jpl/CreateFemModel.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/CreateFemModel.m (revision 3416) +++ /issm/trunk/src/m/solutions/jpl/CreateFemModel.m (revision 3417) @@ -9,8 +9,8 @@ displaystring(md.verbose,'\n reading data from model %s...',md.name); - [m.elements,m.nodes,m.constraints,m.loads,m.materials,m.parameters]=ModelProcessor(md); + [m.elements,m.nodes,m.vertices,m.constraints,m.loads,m.materials,m.parameters]=ModelProcessor(md); displaystring(md.verbose,'%s',' generating degrees of freedom...'); - [m.nodes,m.part,m.tpart]=Dof(m.elements,m.nodes,m.parameters); + [m.nodes,m.vertices,m.part,m.tpart]=Dof(m.elements,m.nodes,m.vertices, m.parameters); displaystring(md.verbose,'%s',' generating single point constraints...'); @@ -30,5 +30,5 @@ displaystring(md.verbose,'%s',' configuring element and loads...'); - [m.elements,m.loads,m.nodes,m.parameters] = ConfigureObjects( m.elements, m.loads, m.nodes, m.materials,m.parameters); + [m.elements,m.loads,m.nodes,m.vertices,m.parameters] = ConfigureObjects( m.elements, m.loads, m.nodes, m.materials,m.parameters); displaystring(md.verbose,'%s',' processing parameters...'); Index: /issm/trunk/src/mex/Dof/Dof.cpp =================================================================== --- /issm/trunk/src/mex/Dof/Dof.cpp (revision 3416) +++ /issm/trunk/src/mex/Dof/Dof.cpp (revision 3417) @@ -12,4 +12,5 @@ /*input datasets: */ DataSet* nodes=NULL; + DataSet* vertices=NULL; DataSet* elements=NULL; DataSet* params=NULL; @@ -28,8 +29,9 @@ FetchData(&elements,ELEMENTS); FetchData(&nodes,NODESIN); + FetchData(&vertices,VERTICESIN); FetchData(¶ms,PARAMS); /*!Generate internal degree of freedom numbers: */ - Dofx(&partition, &tpartition,elements,nodes,params); + Dofx(&partition, &tpartition,elements,nodes, vertices, params); /*partition and tpartition should be incremented by 1: */ @@ -39,4 +41,5 @@ /*write output datasets: */ WriteData(NODES,nodes); + WriteData(VERTICES,vertices); WriteData(PARTITION,partition); WriteData(TPARTITION,tpartition); @@ -44,4 +47,5 @@ /*Free ressources: */ delete nodes; + delete vertices; delete elements; delete params; Index: /issm/trunk/src/mex/Dof/Dof.h =================================================================== --- /issm/trunk/src/mex/Dof/Dof.h (revision 3416) +++ /issm/trunk/src/mex/Dof/Dof.h (revision 3417) @@ -19,16 +19,18 @@ #define ELEMENTS (mxArray*)prhs[0] #define NODESIN (mxArray*)prhs[1] -#define PARAMS (mxArray*)prhs[2] +#define VERTICESIN (mxArray*)prhs[2] +#define PARAMS (mxArray*)prhs[3] /* serial output macros: */ #define NODES (mxArray**)&plhs[0] -#define PARTITION (mxArray**)&plhs[1] -#define TPARTITION (mxArray**)&plhs[2] +#define VERTICES (mxArray**)&plhs[1] +#define PARTITION (mxArray**)&plhs[2] +#define TPARTITION (mxArray**)&plhs[3] /* serial arg counts: */ #undef NLHS -#define NLHS 3 +#define NLHS 4 #undef NRHS -#define NRHS 3 +#define NRHS 4 Index: /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp =================================================================== --- /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp (revision 3416) +++ /issm/trunk/src/mex/ModelProcessor/ModelProcessor.cpp (revision 3417) @@ -14,4 +14,5 @@ DataSet* elements=NULL; DataSet* nodes=NULL; + DataSet* vertices=NULL; DataSet* constraints=NULL; DataSet* loads=NULL; @@ -32,9 +33,10 @@ /*Create elements, nodes and materials: */ - CreateDataSets(&elements,&nodes,&materials,&constraints, &loads, ¶meters, iomodel,MODEL); + CreateDataSets(&elements,&nodes,&vertices,&materials,&constraints, &loads, ¶meters, iomodel,MODEL); /*Write output data: */ WriteData(ELEMENTS,elements); WriteData(NODES,nodes); + WriteData(VERTICES,vertices); WriteData(CONSTRAINTS,constraints); WriteData(LOADS,loads); @@ -47,4 +49,5 @@ delete elements; delete nodes; + delete vertices; delete constraints; delete loads; Index: /issm/trunk/src/mex/ModelProcessor/ModelProcessor.h =================================================================== --- /issm/trunk/src/mex/ModelProcessor/ModelProcessor.h (revision 3416) +++ /issm/trunk/src/mex/ModelProcessor/ModelProcessor.h (revision 3417) @@ -26,12 +26,13 @@ #define ELEMENTS (mxArray**)&plhs[0] #define NODES (mxArray**)&plhs[1] -#define CONSTRAINTS (mxArray**)&plhs[2] -#define LOADS (mxArray**)&plhs[3] -#define MATERIALS (mxArray**)&plhs[4] -#define PARAMETERS (mxArray**)&plhs[5] +#define VERTICES (mxArray**)&plhs[2] +#define CONSTRAINTS (mxArray**)&plhs[3] +#define LOADS (mxArray**)&plhs[4] +#define MATERIALS (mxArray**)&plhs[5] +#define PARAMETERS (mxArray**)&plhs[6] /* serial arg counts: */ #undef NLHS -#define NLHS 6 +#define NLHS 7 #undef NRHS #define NRHS 1