Index: /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 3424) +++ /issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 3425) @@ -2,5 +2,4 @@ * CreateElementsNodesAndMaterialsSlopeCompute.c: */ - #include "../../DataSet/DataSet.h" @@ -15,12 +14,7 @@ void CreateElementsNodesAndMaterialsSlopeCompute(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: */ DataSet* elements = NULL; @@ -29,66 +23,4 @@ DataSet* materials = NULL; - /*Objects: */ - Node* node = NULL; - Vertex* vertex = NULL; - Tria* tria = NULL; - Penta* penta = 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 - - /*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_s[3]; - double tria_b[3]; - bool tria_onwater; - - /*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_s[6]; - double penta_b[6]; - int penta_onbed; - bool penta_onwater; - - /* 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()); @@ -97,41 +29,6 @@ materials = new DataSet(MaterialsEnum()); - /*Now, is the flag isstokes on? otherwise, do nothing: */ - if (!iomodel->isstokes & !iomodel->ishutter)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); - - - /*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: */ @@ -146,58 +43,12 @@ for (i=0;inumberofelements;i++){ - #ifdef _PARALLEL_ - /*!All elements have been partitioned above, only create elements for this CPU: */ - if(my_rank==epart[i]){ - #endif - - - /*ids: */ - tria_id=i+1; //matlab indexing. - tria_matice_id=-1; //no materials - tria_matpar_id=-1; //no materials - tria_numpar_id=1; //no materials + if(iomodel->my_elements[i]){ + + /*Create and add tria element to elements dataset: */ + elements->AddObject(new Tria(i,iomodel)); - /*vertices offsets: */ - 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); - - /*surface and bed:*/ - 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)); - - /*element on water? : */ - tria_onwater=(bool)*(iomodel->elementonwater+i); - - /*Create properties: */ - tria_properties=new ElementProperties(3,NULL, tria_s, tria_b, NULL, NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, UNDEF, 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); - - #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 + /*Create and add material property to materials dataset: */ + materials->AddObject(new Matice(i,iomodel,3)); + } }//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]){ + /*Create and add penta element to elements dataset: */ + elements->AddObject(new Penta(i,iomodel)); - - /*name and id: */ - penta_id=i+1; //matlab indexing. - penta_matice_id=-1; //no materials - penta_matpar_id=-1; //no materials - penta_numpar_id=1; //no materials - + /*Create and add material property to materials dataset: */ + materials->AddObject(new Matice(i,iomodel,6)); - /*vertices,thickness,surface,bed and drag: */ - for(j=0;j<6;j++){ - penta_node_ids[j]=(int)*(iomodel->elements+elements_width*i+j); - 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)); - } - - /*diverse: */ - penta_onbed=(int)*(iomodel->elementonbed+i); - penta_onwater=(bool)*(iomodel->elementonwater+i); - - /*Create properties: */ - penta_properties=new ElementProperties(6,NULL, penta_s, penta_b, NULL, NULL, NULL, NULL, UNDEF, UNDEF, UNDEF,UNDEF,penta_onbed, penta_onwater, UNDEF, UNDEF, 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); - - #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;inumberofnodes); - - 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){ + /*Add new constrant material property tgo materials, at the end: */ + materials->AddObject(new Matpar(iomodel)); - /*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){ @@ -364,78 +109,17 @@ IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); + for (i=0;inumberofvertices;i++){ - /*Get number of dofs per node: */ - DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); + /*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)); - for (i=0;inumberofnodes;i++){ - #ifdef _PARALLEL_ - /*keep only this partition's nodes:*/ - if((my_grids[i]==1)){ - #endif + /*Add node to nodes dataset: */ + nodes->AddObject(new Node(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); - - - 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); - - /*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->gridonbed[i]==0){ - 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(); - vertices->Presort(); - materials->Presort(); /*Clean fetched data: */ @@ -452,16 +136,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: