Index: /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp (revision 3427) +++ /issm/trunk/src/c/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp (revision 3428) @@ -14,11 +14,6 @@ void CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ - /*output: int* epart, int* my_grids, double* my_bordergrids*/ - - int i,j,k,n; - extern int my_rank; - extern int num_procs; /*DataSets: */ @@ -28,85 +23,4 @@ DataSet* materials = NULL; - /*Objects: */ - Node* node = NULL; - Vertex* vertex = NULL; - Penta* penta = NULL; - 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: */ - int elements_width; //size of elements - double B_avg; - - /*matice constructor input: */ - int matice_mid; - double matice_B; - double matice_n; - - /*penta constructor input: */ - - int penta_id; - int penta_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_geothermalflux[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 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; - /*First create the elements, nodes and material properties: */ elements = new DataSet(ElementsEnum()); @@ -115,25 +29,6 @@ materials = new DataSet(MaterialsEnum()); - /*return if 2d mesh*/ - if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return; - - /*Width of elements: */ - elements_width=6; //penta elements - - #ifdef _PARALLEL_ - /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ - 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->elements2d); - - /*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); /*Fetch data needed: */ @@ -155,92 +50,12 @@ 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->my_elements[i]){ - - /*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; + /*Create and add tria element to elements dataset: */ + elements->AddObject(new Tria(i,iomodel)); - /*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_geothermalflux[j]=*(iomodel->geothermalflux+ ((int)*(iomodel->elements+elements_width*i+j)-1)); + /*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); - - /*We need the field collapse for transient, so that we can use compute B with the average temperature*/ - 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 properties: */ - penta_properties=new ElementProperties(6,penta_h, penta_s, penta_b, penta_k, NULL, NULL, penta_geothermalflux, 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 - }//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 - - /*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: */ + /*Create nodes and vertices: */ if (strcmp(iomodel->meshtype,"3d")==0){ IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); @@ -362,69 +95,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]){ + + /*Add vertex to vertices dataset: */ + vertices->AddObject(new Vertex(i,iomodel)); - /*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); + /*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 (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,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: */ @@ -441,4 +122,11 @@ xfree((void**)&iomodel->gridoniceshelf); + /*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: @@ -449,15 +137,3 @@ *pmaterials=materials; - /*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 - }