Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3426) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3427) @@ -15,112 +15,12 @@ void CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ + /*Intermediary*/ int i,j,k,n; - extern int my_rank; - extern int num_procs; /*DataSets: */ DataSet* elements = NULL; DataSet* nodes = NULL; - DataSet* vertices = 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; - - /*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]; - 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_mid; - int penta_mparid; - int penta_numparid; - int penta_g[6]; - double penta_h[6]; - double penta_s[6]; - double penta_b[6]; - double penta_k[6]; - 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]; - int penta_thermal_steadystate; - 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 vertex_partitionborder=0; - 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; - /*Nodes cloning*/ - double e1,e2; - int i1,i2; - int pos; - #endif - int first_grid_index; /*First create the elements, nodes and material properties: */ @@ -130,37 +30,8 @@ materials = new DataSet(MaterialsEnum()); - /*Width of elements: */ - if(strcmp(iomodel->meshtype,"2d")==0){ - elements_width=3; //tria elements - } - else{ - ISSMERROR("not implemented yet!"); - 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(3*iomodel->numberofelements,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); /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ - /*2d mesh: */ if (strcmp(iomodel->meshtype,"2d")==0){ @@ -176,63 +47,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 need for materials - tria_matpar_id=-1; //no need for materials - tria_numpar_id=1; + if(iomodel->my_elements[i]){ - /*vertices offsets: */ - tria_node_ids[0]=3*i+1; - tria_node_ids[1]=3*i+2; - tria_node_ids[2]=3*i+3; + /*Create and add tria element to elements dataset: */ + elements->AddObject(new Tria(i,iomodel)); - /*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)); - - /*element on iceshelf?:*/ - 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, NULL, NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, 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); - - #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[3*i+0]=1; - my_grids[3*i+1]=1; - my_grids[3*i+2]=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;iAddObject(new Matpar(iomodel)); - /*Get edges and elements*/ - IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges"); - - /*!All elements have been partitioned above, only create elements for this CPU: */ - for (i=0;inumberofedges;i++){ - - /*Get left and right elements*/ - e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2] - e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2] - - /* 1) If the element e1 is in the current partition - * 2) and if the edge of the element is shared by another element - * 3) and if this element is not in the same partition: - * we must clone the nodes on this partition so that the loads (Numericalflux) - * will have access to their properties (dofs,...)*/ - if(my_rank==epart[(int)e1] && !isnan(e2) && my_rank!=epart[(int)e2]){ - - /*1: Get vertices ids*/ - i1=(int)iomodel->edges[4*i+0]; - i2=(int)iomodel->edges[4*i+1]; - - /*2: Get the column where these ids are located in the index*/ - pos==UNDEF; - for(j=0;j<3;j++){ - if (iomodel->elements[3*(int)e2+j]==i1) pos=j+1; - } - ISSMASSERT(pos!=UNDEF); - - /*3: We have the id of the elements and the position of the vertices in the index - * we can now create the corresponding nodes:*/ - my_grids[(int)e2*3+pos-1]=1; - my_grids[(int)e2*3+((pos+1)%3)]=1; - } - } - - /*Free data: */ - xfree((void**)&iomodel->edges); - - /*From the element partitioning, we can determine which nodes are on the inside of this cpu's - *element partition, and which are on its border with other nodes:*/ - gridborder=NewVec(3*iomodel->numberofelements); - - for (i=0;i<3*iomodel->numberofelements;i++){ - if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); - } - VecAssemblyBegin(gridborder); - VecAssemblyEnd(gridborder); - - VecToMPISerial(&my_bordergrids,gridborder); - - #endif - - /*Add one constant material property to materials: */ - matpar_mid=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); - - /*Ok, let's summarise. Now, every CPU has the following array: my_grids - 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"); - } + /*Create nodes and vertices: */ IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x"); IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y"); @@ -346,83 +80,54 @@ IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed"); IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface"); + IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter"); 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); - - /*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 + if (strcmp(iomodel->meshtype,"3d")==0){ + IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); + IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); } - /*Build Nodes dataset -> 3 for each element*/ + /*Build Vertices dataset*/ + for (i=0;inumberofvertices;i++){ + + /*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)); + + } + } + + /*Build Nodes dataset -> 3 for each element (Discontinuous Galerkin)*/ for (i=0;inumberofelements;i++){ for (j=0;j<3;j++){ - #ifdef _PARALLEL_ - /*!All elements have been partitioned above, only create elements for this CPU: */ - if(my_grids[3*i+j]){ - #endif + if(my_nodes[3*i+j]){ //Get id of the vertex on which the current node is located k=(int)*(iomodel->elements+elements_width*i+j)-1; //(Matlab to C indexing) - ISSMASSERT(k>=0 && knumberofnodes); - - //Get node properties - node_id=i*3+j+1; - #ifdef _PARALLEL_ - if(my_bordergrids[node_id-1]>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[k], - (int)iomodel->gridonsurface[k], - (int)iomodel->gridoniceshelf[k], - (int)iomodel->gridonicesheet[k]); - - if (strcmp(iomodel->meshtype,"3d")==0){ - if (isnan(iomodel->uppernodes[k])){ - node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves. - } - else{ - node_upper_node_id=(int)iomodel->uppernodes[k]; - } - } - 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, k, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties); + ISSMASSERT(k>=0 && knumberofvertices); /*Add node to nodes dataset: */ - nodes->AddObject(node); - #ifdef _PARALLEL_ + nodes->AddObject(new Node(k,i,iomodel)); + } - #endif } } + + /*Clean fetched data: */ + xfree((void**)&iomodel->deadgrids); + xfree((void**)&iomodel->x); + xfree((void**)&iomodel->y); + xfree((void**)&iomodel->z); + xfree((void**)&iomodel->thickness); + xfree((void**)&iomodel->bed); + xfree((void**)&iomodel->gridonbed); + xfree((void**)&iomodel->gridonsurface); + xfree((void**)&iomodel->gridonhutter); + xfree((void**)&iomodel->uppernodes); + xfree((void**)&iomodel->gridonicesheet); + 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 @@ -432,30 +137,4 @@ vertices->Presort(); materials->Presort(); - - /*Clean fetched data: */ - xfree((void**)&iomodel->deadgrids); - xfree((void**)&iomodel->elements); - xfree((void**)&iomodel->x); - xfree((void**)&iomodel->y); - xfree((void**)&iomodel->z); - xfree((void**)&iomodel->thickness); - xfree((void**)&iomodel->bed); - xfree((void**)&iomodel->gridonbed); - xfree((void**)&iomodel->gridonsurface); - xfree((void**)&iomodel->uppernodes); - 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: Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 3426) +++ /issm/trunk/src/c/objects/Node.cpp (revision 3427) @@ -52,5 +52,5 @@ } /*}}}*/ -/*FUNCTION Node constructor from iomodel{{{2*/ +/*FUNCTION Node constructor from iomodel continuous Galerkin{{{2*/ Node::Node(int i, IoModel* iomodel){ //i is the node index @@ -131,4 +131,51 @@ } } +} +/*}}}*/ +/*FUNCTION Node constructor from iomodel discontinuous Galerkin{{{2*/ +Node::Node(int i,int j,IoModel* iomodel){ + /* i -> index of the vertex in C indexing + * j -> index of the node in C indexing*/ + + int numdofs; + int partitionborder; + int vertex_id; + int upper_node_id; + + /*id: */ + this->id=j+1; //matlab indexing + + /*indexing:*/ + DistributeNumDofs(&numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); //number of dofs per node + if(iomodel->my_bordervertices[i])partitionborder=1; else partitionborder=0;//is this node on a partition border? + + this->indexing.Init(numdofs,partitionborder); + + /*properties (come from vertex number i): */ + this->properties.Init( + (int)iomodel->gridonbed[i], + (int)iomodel->gridonsurface[i], + (int)iomodel->gridoniceshelf[i], + (int)iomodel->gridonicesheet[i]); + + /*hooks: */ + vertex_id=i+1; //matlab indexing + + 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); + this->hupper_node.Init(&upper_node_id,1); + } /*}}}*/ Index: /issm/trunk/src/c/objects/Node.h =================================================================== --- /issm/trunk/src/c/objects/Node.h (revision 3426) +++ /issm/trunk/src/c/objects/Node.h (revision 3427) @@ -44,4 +44,5 @@ Node(int id,DofIndexing* indexing, NodeProperties* properties, Hook* vertex, Hook* upper_node); Node(int i, IoModel* iomodel); + Node(int i,int j,IoModel* iomodel); ~Node(); /*}}}*/