Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 3420) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 3421) @@ -12,8 +12,7 @@ #include "../IoModel.h" - void CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ - + /*Intermediary*/ int i,j,k,n; @@ -146,5 +145,4 @@ } //if (strcmp(meshtype,"2d")==0) - /*Add new constrant material property tgo materials, at the end: */ @@ -166,5 +164,4 @@ IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); } - for (i=0;inumberofvertices;i++){ Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 3420) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 3421) @@ -2,5 +2,4 @@ * CreateElementsNodesAndMaterialsDiagnosticVert.c: */ - #include "../../DataSet/DataSet.h" @@ -13,14 +12,8 @@ #include "../IoModel.h" - void CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ - - /*output: int* epart, int* my_grids, double* my_bordergrids*/ - - + /*Intermediary*/ int i,j,k,n; - extern int my_rank; - extern int num_procs; /*DataSets: */ @@ -36,88 +29,7 @@ Matice* matice = NULL; Matpar* matpar = 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: */ - const int elements_width=6; - double B_avg; - - /*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]; - int penta_shelf; - int penta_onbed; - int penta_onsurface; - int penta_collapse; - double penta_melting[6]; - double penta_accumulation[6]; - int penta_artdiff; - double penta_viscosity_overshoot; - double penta_stokesreconditioning; - 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 vertex_id; - int node_partitionborder=0; - int node_onbed; - int node_onsurface; - int node_onshelf; - int node_onsheet; - 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; - - /*Penalty partitioning: */ - int num_grids3d_collapsed; - double* double_penalties_grids3d_collapsed=NULL; - double* double_penalties_grids3d_noncollapsed=NULL; - int grid_ids[6]; - int num_grid_ids; - int grid_id; - - + /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ + if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return; /*First create the elements, nodes and material properties: */ @@ -127,28 +39,6 @@ materials = new DataSet(MaterialsEnum()); - /*Now, is the iomodel running in 3d? : */ - if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return; - - #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); - - /*Used later on: */ - my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int)); - #endif - + /*Partition elements and vertices and nodes: */ + Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle); /*Create 3d elements: */ @@ -167,71 +57,14 @@ for (i=0;inumberofelements;i++){ - #ifdef _PARALLEL_ - /*We are using our element partition to decide which elements will be created on this node: */ - if(my_rank==epart[i]){ - #endif - - /*name and id: */ - penta_id=i+1; //matlab indexing. - penta_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; + if(iomodel->my_elements[i]){ - /*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_melting[j]=*(iomodel->melting+ ((int)*(iomodel->elements+elements_width*i+j)-1))/iomodel->yts; - penta_accumulation[j]=*(iomodel->accumulation+ ((int)*(iomodel->elements+elements_width*i+j)-1))/iomodel->yts; + /*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)); + } - - /*diverse: */ - penta_shelf=(int)*(iomodel->elementoniceshelf+i); - penta_onbed=(int)*(iomodel->elementonbed+i); - penta_onsurface=(int)*(iomodel->elementonsurface+i); - penta_collapse=1; - penta_onwater=(bool)*(iomodel->elementonwater+i); - - /*Create properties: */ - penta_properties=new ElementProperties(6,penta_h, penta_s, penta_b, NULL, penta_melting, penta_accumulation, NULL, UNDEF,UNDEF,UNDEF, 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. - - /*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 }//for (i=0;ielementonwater); - - /*Add one constant material property to materials: */ - 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; - - /*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); + /*Add new constrant material property tgo materials, at the end: */ + materials->AddObject(new Matpar(iomodel)); - #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 - - - /*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"); - } + IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); + IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x"); IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y"); @@ -321,64 +99,17 @@ 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;inumberofvertices;i++){ - for (i=0;inumberofnodes;i++){ - #ifdef _PARALLEL_ - /*keep only this partition's nodes:*/ - if((my_grids[i]==1)){ - #endif + /*vertices and nodes (same number, as we are running continuous galerkin formulation: */ + if(iomodel->my_vertices[i]){ - /*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])); + /*Add vertex to vertices dataset: */ + vertices->AddObject(new Vertex(i,iomodel)); - vertices->AddObject(vertex); + /*Add node to nodes dataset: */ + nodes->AddObject(new Node(i,iomodel)); - - 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_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. - } - else{ - node_upper_node_id=(int)iomodel->uppernodes[i]; - } - - /*Create node using its constructor: */ - node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties); - - /*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(); - vertices->Presort(); - materials->Presort(); /*Clean fetched data: */ @@ -395,16 +126,10 @@ 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 + /*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(); cleanup_and_return: